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INTRODUCTION 


The  work  under  this  grant  on  plasma  presheaths,  which  form  a  transition  region  be¬ 
tween  the  coliisionless  electrode  sheaths  and  the  plasma,  is  directed  toward  the  problems 
of  the  Thermionic  Energy  Convertor  (TEC).  Figure  1  shows  a  schematic  of  a  TEC  in  a 
reactor  core  for  space  power  applications  and  the  basic  physics.  Cesium  is  put  the  gap 
between  the  emitter  and  collector  for  two  purposes:  first,  to  ionize  and  neutralize  the 
space  charge  so  that  a  useful  electron  current  density  can  pass  (10  -  100  amps/square 
cm),  and  second  to  reduce  the  electrode  work  functions  by  adsorption  of  cesium.  Of  the 
plasma  physics  of  the  the  cesium  filled  gap  of  the  TEC,  the  plasma-electrode  interactions 
are  the  most  significant  part  because  these  regions  form  boundary  conditions  which  con¬ 
trol  the  plasma  density  and  temperatures  of  the  entire  gap.  Thus  the  research  under  this 
grant  has  been  directed  toward  the  study  of  collisional  presheaths  which  form  the  layer 
adjacent  to  an  electrode  on  the  order  of  one  ion  mean  free  path  thick.  However,  the  re¬ 
search  pursued  under  this  grant  is  not  limited  in  applicabilty  to  TECs  but  is  of  interest 
to  plasma-surface  interactions  in  general.  Other  applications  include  electric  propulsion 
where  electrode  erosion  is  a  problem  and  not  fully  understood  and  more  generally  any 
plasmarsurface  interaction. 

This  report  includes  the  asymptotic  presheath  theory  developed,  and  is  preceded  by  the 
basic  theory  of  the  Thermionic  Energy  Convertor  (TEC)  and  is  followed  by  the  application 
of  the  theory  to  a  time  dependent  model  of  the  TEC  in  the  program  called  TEC.  As  shown 
in  the  TEC  results,  the  agreement  with  experiment  is  good  except  in  the  low  current  regime 
of  the  TEC  where  an  unexplained  disagr*-  w.t  remains.  This  is  still  a  puzzie. 
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BASIC  TEC  THEORY 


The  basic  theory  of  the  TEC  is  set  forth  in  the  following  paper  published  under  this 
grant. 


4 


T-PS/ 15/3/13329 


Effects  of  Emitter  Sheath  Ion  Reflection  and  Trapped 
Ions  on  Thermionic  Converter  Performance 
Using  an  Isothermal  Electron  Model 

G.  L.  Main 
S.  H.  Lam 


Reprised  from:  IEEE  TRANSACTIONS  ON  PLASMA  SCIENCE.  Vol.  PS-13,  No.  3.  JUNE  1917. 


UK  TRANSACTIONS  ON  PLASMA  SCIENCE.  VOL.  PS-13,  NO.  ).  JUNE  I9S7 


Effects  of  Emitter  Sheath  Ion  Reflection  and  T  rapped 
Ions  on  Thermionic  Converter  Performance  Using  an 

Isothermal  Electron  Model 

GEOFFREY  L.  MAIN  and  S.  H.  LAM 


yU«Mcf— -TMi  paper  cowptos  exact  coBblaakw  ffecatb  caJcalatJoa* 
ta  aa  bMbtn—t  ctoctroq  nSrl  of  a  tkwlaak  caavrrtrr.  The  vtnlttvr 
sbeatb  Nnactar*  takas  lato  account  reflected  leas,  trapped  too*.  aad 
seriate  eiateilea  iaas.  It  is  shows  that  tesetaiag  the  act  loss  of  loos  at 
the  emitter  la  the  ifaited  saode  by  these  pheasaieaa  defrades  perfor- 
auace.  la  addkioa,  h  b  shows  that  whea  the  eadtter  returns  too  many 
of  the  iaas,  the  ate  b  eating  ebbed  because  there  b  Insufficient  resistive 
heating  to  eiaiaUia  the  necessary  plasma  electron  temperature  for  ion* 
batba.  These  resnlu  suggest  that  the  ignited  mode  cannot  be  improved 
much.  However,  aoaignited  modes  in  which  the  electron  temperature 
remains  tow,  suck  as  the  poised  mode,  do  aat  safer  from  thb  adverse 
behavior. 


I.  Introduction 

EMITTER  sheath  phenomena  are  important  in  ther¬ 
mionic  energy  converters  because  the  emitter  sheath 
forms  the  emitter  boundary  condition  for  the  plasma  in 
the  gap  by  controlling  both  the  ion  loss  rate  and  the  loss 
rate  of  hot  ( 3000  K )  plasma  electrons  to  the  emitter.  This 
paper  examines  two  expected  emitter  sheath  phenomena 
and  their  effects  on  converter  performance:  reflection  of 
ions  coming  from  the  plasma  by  a  double  emitter  sheath, 
and  ions  trapped  in  the  double  emitter  sheath.  The  authors 
have  previously  suggested  that  ion  reflection  might  im¬ 
prove  thermionic  energy  converter  performance  ( 1]  and 
have  subsequently  shown  that  ion  reflection  at  the  emitter 
is  likely  to  degrade  the  performance  in  the  ignited  mode 
and,  in  addition,  that  trapped  ions  in  a  double  emitter 
sheath  are  also  likely  to  degrade  performance  in  the  ig¬ 
nited  mode  (2).  Lundgren  [3],  [4]  has  also  shown  this  with 
simplified  km  and  electron  dynamics.  In  the  present  paper 
the  effects  of  emitter  ion  reflection  and  ion  trapping  in  the 
ignited  mode  are  calculated  using  exact  electron  and  ion 
dynamics  in  the  collisionless  (except  for  ion  trapping) 
sheaths.  The  electrons  entering  the  sheaths  from  the 
plasma  are  assumed  to  have  a  Maxwellian  distribution, 
but  no  assumptions  are  made  about  the  returning  elec¬ 
tions,  and  the  electron  density  in  the  sheath  is  calculated 
exactly.  The  ions  entering  the  sheaths  from  the  plasma  are 
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not  assumed  cold,  but  are  given  the  correct  ion  tempera¬ 
ture  and  shifted  in  velocity  according  io  a  generalization 
of  the  Bohm  criterion  (5],  {6J. 

Both  ion  reflection  and  trapped  ions  in  the  emitter  sheath 
reduce  the  normalized  (by  plasma  density)  net  ion  loss 
rate  to  the  emitter.  Also,  both  of  these  phenomena  raise 
the  normalized  plasma  density  adjacent  to  the  emitter.  The 
higher  plasma  density  at  the  emitter  causes  a  greater  in¬ 
crease  in  the  loss  of  hot  plasma  electron  energy  to  the 
emitter  than  the  corresponding  decrease  in  the  loss  of  ion¬ 
ization  energy  (carried  by  the  ions)  to  the  emitter.  There¬ 
fore,  these  emitter  sheath  phenomena  increase  arc-drop 
Within  the  limitations  of  the  present  isothermal  ther¬ 
mionic  converter  formulation,  all  three  of  these  phenom¬ 
ena  (which  become  significant  at  low  currents)  steepen  the 
current-voltage  characteristic.  At  low  current  densities, 
the  present  theory  shows  that  the  collector  sheath  height 
decreases,  resulting  in  a  larger  electron  diffusion  velocity 
than  can  be  justified  for  the  continuum  model  used  in  the 
plasma  region.  The  result  of  lower  performance  at  lower 
current  is  in  agreement  with  experimental  studies.  At  some 
current  density  which  depends  strongly  on  the  emitter 
sheath  conditions,  the  ignited  mode  is  no  longer  self-sus¬ 
taining  and  the  arc  is  extinguished. 

Fig-  1  is  a  schematic  diagram  of  the  cesium  diode  con¬ 
verter.  The  emitter  is  heated  externally  to  temperature  7> 
which  is  typically  1750  K  or  higher,  and  the  collector  is 
cooled  to  temperature  Tc  which  is  typically  900-1 100  K. 
The  gap  space  d,  or  convener  length,  which  is  typically 
0.25  mm,  separates  the  emitter  from  the  collector.  The 
cesium  reservoir,  which  is  sometimes  imbedded  in  the 
collector,  is  kept  at  temperature  7*  to  maintain  the  desired 
cesium  pressure  (typically  1  to  2  ton)  in  the  gap.  The 
electrical  load  is  connected  across  the  emitter  and  collec¬ 
tor  to  produce  power. 

II.  The  Isothermal  Electron  Formulation 

In  this  section  the  isothermal  thermionic  convener  for¬ 
mulation  is  developed.  The  formulation  is  similar  to  that 
of  Lam  [7}  but  is  generalized  to  eliminate  the  assumption 
of  high  sheaths  which  has  previously  been  used  to  sim¬ 
plify  the  electron  dynamics.  Since  both  low-emitter  and 
low-collector  sheath  heights  are  encountered  as  a  conse¬ 
quence  of  ion  reflection  and  trapped  ions,  the  assumption 
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Fijt  I .  The  cesium  diode  convener. 


«f  Boltzmann  plasma  electron  distributions  at  the  plasma- 
sheath  interface  must  be  abandoned.  At  both  the  emitter 
and  collector,  the  low  sheaths  return  few  plasma  elec- 
tnms.  leaving  the  distributions  largely  one  sided.  Fur¬ 
thermore.  at  the  emitter  sheath  emitted  electrons  must  be 
taken  into  account.  Thus  the  ratio  of  electrons  moving 
toward  the  sheath  to  the  total  density  of  electrons  at  the 
sheath  edge  is  not  I  /2.  as  in  the  Boltzmann  assumption. 

In  Fig.  2  we  define  the  potentials  in  the  converter.  All 
of  the  potentials  are  nondimensionalized  by  emitter  tem¬ 
perature  as  follows: 


* 


q4> 

kTt 


(0 


where 


The  emitted  current  density  which  crosses  the  emitter 
sheath  potential  peak  into  the  converter  plasma  region  is 

Jt  -  7*  exp  (-Ax).  Ax  >  0 

J£-J„  Ax  SO.  (5) 

We  also  define  the  net  current  density  through  the  con¬ 
verter  J  and  the  normalized  current  density 


We  have  assumed  for  convenience  that  the  ion  contribu¬ 
tion  to  net  current  is  negligible  because  the  cesium-ion- 
to-electron-mass  ratio  is  enormous.  Ions  will  typically 
contribute  no  more  than  I  percent  of  the  net  current.  Elec¬ 
tron  temperature  is  nondimensionalized  as 


X  nondimensional  potential. 

6  potential. 

*1  electron  charge, 
k  Boltzmann  constant,  and 
7V  emitter  temperature. 

Wc  also  use  the  following  terminology  for  various  poten¬ 
tials  in  the  converter: 

♦»  emitter  work  function. 

Ay  back  sheath  height. 

A\,  reflective  potential. 

Xj  emitter  sheath  height. 

A\r  plasma  potential  drop. 

Vj  arc-drop. 

\,  collector  sheath  height. 

♦,  collector  work  function,  and 

l  .„,  converter  output  voltage. 

Inspection  of  Fig.  2  immediately  yields  the  following  re¬ 
lations: 

*0  -  AX  (2) 

VJ  “  (Xc  “  X*)  ~  A xr.  (3) 

The  Richardson  current  density  of  electrons  from  the 
emitter  is 

A  (A)  -  120  n(K:)  exp  (-♦*).  (4) 


where  T,  is  the  plasma  electron  temperature  which,  in  this 
section,  is  constant  by  the  isothermal  assumption.  Fi¬ 
nally,  we  have  the  thermal  speeds: 


The  isothermal  formulation  is  developed  from  here  in 
the  same  way  as  the  general  formulation  except  that  we 
take  full  advantage  of  the  isothermal  assumption  by  look¬ 
ing  only  at  the  global  conservation  equations  instead  of 
the  local  ones  used  in  the  general  formulation.  We  then 
assume  that  the  transport  properties,  collision  frequen¬ 
cies.  and  the  ionization  source  coefficient  are  constant 
across  the  convener  because  of  the  isothermal  assump¬ 
tion.  Also,  we  find  only  the  steady-state  solution.  We 
carry  out  this  development  by  deriving  the  global  conser¬ 
vation  equations  for  the  isothermal  case  (current,  momen¬ 
tum.  and  electron  energy)  and  then  reducing  these  to  a  set 
of  three  simultaneous  equations  in  the  variables  r,  xe<  and 
Xc-  In  some  cases  the  actual  calculations  are  carried  out 
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using  different  variables  when  \t  or  Xc  are  small  or  zero. 
!n  the  case,  for  instance,  of  a  single  ion-repelling  emitter 
sheath  we  use  j  because  Xf.  is  zero.  These  equations  are 
nonlinear  and  solved  numerically  using  a  positive  definite 
Newton's  method. 

First,  we  consider  conservation  of  current.  The  collec¬ 
tor  is  assumed  to  emit  nothing:  therefore,  at  the  plasma- 
collector  sheath  interface  we  have 


J  = 


a,.a\n(  I ) 

2 


/» 


(10) 


where  a,  is  the  fraction  of  the  total  plasma  density  at  the 
collector  sheath  which  is  moving  toward  the  collector  and 
n  ( 1 )  is  the  total  plasma  density  at  the  plasma-collector 
sheath  interface.  Because  we  continue  to  assume  that  the 
pan  of  the  plasma  electron  distribution  coming  into  the 
collector  sheath  is  Maxwellian,  we  can  write  a,  as 


1 


2  f"*'T  - 

'  +  Tr  J.  dU 


(ID 


which  takes  into  account  the  plasma  electrons  reflected  by 
the  collector  sheath.  We  still  assume  that  the  plasma  elec¬ 
tron  distribution  coming  into  the  collector  sheath  is  Max¬ 
wellian  and  that  it  docs  not  have  any  velocity  shift  be¬ 
cause  the  sheath  is  expected  to  be  electron  repelling.  In 
the  limit  of  a  high  collector  sheath,  a,  *  1/2  and  we 
have  a  fully  Boltzmann  distribution  of  electrons  at  the  col¬ 
lector  sheath  edge.  The  situation  at  the  emitter  is  more 
complex  because  the  emitted  electrons  must  be  taken  into 
account.  We  have  the  backscattcred  current  density  JBS 
which  is  the  plasma  electron  current  density  moving  into 
the  emitter: 


where 


is  the  electron  Mach  number  at  the  emitter.  This  is  |ust 
an  application  of  ( 13). 

Electron  energy  conservation  is  developed  by  consid¬ 
ering  energy  exchange  with  the  emitter  and  collector  and 
energy  lost  to  ionization.  Power  carried  into  the  plasma 
by  emitted  electrons  is 

Ph  =  Jt  ( 2  +  +  Ax)  —  (17) 

Power  returned  to  the  emitter  is 

Pbs-  Uf  -J)(2t  +  ♦,  +  Ax)— .  (18) 

</ 

Power  flowing  into  the  collector  is 

Pc  -  y(2r  +  ♦, .  +  Vj  +  Ax)—.  (19) 

Ionization  power  loss  is 

Pu«  «  -U  ^  —  (20) 

where  JHm  is  the  total  ion  current  into  both  the  emitter  and 
collector,  and  Vp  is  the  first  ionization  energy.  Conser¬ 
vation  of  electron  energy  is 

Pf  -  +  Pc  +  P.*  (21) 

which  can  be  reduced  to 


Jbs 


n(0)grq„  (  Xt:\ 
2  eXP\  r) 


(12) 


where  ir(0)  is  the  total  plasma  density  at  the  emitter 
sheath-plasma  interface  and  a„  is  the  fraction  of  total 
plasma  density  at  the  interlace  moving  toward  the  emitter. 
Continuity  of  electron  current  demands 

h  «  J,s  +  J  (13) 

which  can  be  written  as 


h  «  J  I  + 


n(O)ai) 

«(l)0| 


exp 


(14) 


This  can  be  rewritten  using  (3)  and  (6)  as 

I 


'  .  .  «(0)<V  (V*  +  AxP\ 

1  STTTST exp  V  7  / 

The  quamiry  a,,  can  be  written  as 

-■iHj-'Mr) 


(15) 


(16) 


t»I  -  \jVj  -  \j,Vri  (22) 

where  j,  =  In  the  ignited  mode  r  is  generally 

about  2(77  *  1750  K  and  Tr  -  3000  K).  consequently 
the  arc-drop  K(,  is  negative.  In  other  words,  the  high 
plasma  electron  temperature  is  generated  by  resistance 
heating. 

Finally,  we  consider  electron  and  ion  momentum.  From 
electron  momentum  conservation .  wc  find  the  potential 
drop  in  the  plasma  region.  By  adding  the  electron  and  ion 
momentum  equations  as  in  the  general  case,  wc  find  our 
diffusion  equation  and  boundary  conditions  to  which  the 
sheaths  contribute  flux  terms.  When  we  introduce  the  ion¬ 
ization  source  term  into  this,  wc  have  the  complete  for¬ 
mulation.  Electron  momentum  conservation  is 


0.  MS.  (23) 

dx  dx  \ 

where  X,  is  electron  mean-free  path.  Using  pr  *  nkTr  and 
J  »  qnu„  we  can  rearrange  (23)  into 


J  « 


dn 

7x  +  nq 


dA\ 

dx)' 


t 


«u 


(24) 
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This  can  be  further  reduced  by  dividing  by  Jt  and  using 
|  *  x  / d  where  J  is  the  convener  gap  thickness. 

»  K  I  (  tin  d\  \ 

'  1251 

Integration  of  this  equation  from  the  emitter  sheath  inter¬ 
face  us  the  coliector  sheath  interface  yields 


where 


where 


ClSi)  * J* 

4  i /  r  f'  nt 

r  -  -  r  ^  -77:  ‘it. 


The  quantity  R  is  the  normalized  plasma  resistance. 

The  ion  and  electron  momentum  equations  can  be  writ¬ 
ten 

tin  4 H  mint, a,  (,Q  4 

tT—  «  -qn- - - -  (28a) 

dx  dx  X, 


, ^  tin  dt  Mini, it, 

"■ir*"*-  —  ,28b| 

when:  \  is  ion  mean-free  path  and  a,  is  ion  thermal  speed: 


Bit 

Addition  of  (28a)  and  (28b)  yields 


(if.  +  kTt )  y 
dx 


which  is  amhipolar  diffusion.  Equation  (29)  is  differi-n- 
tiated  to  become 

d'n  a,  m  <1 

(*r.  +  47V)  J--  +  -7“  7  (mm.  ) 

dx '  A,  d  \ 

+  •““(«",)*()•  (30) 

A,  dx 

We  assume  recombination  is  negligible  and  the  ionization 
source  term  is. 

-7  (*.»)  *  -7  (m,m)  *  Sn.  (31) 

dx  dx 


Using  (311  in  (30)  yields 
d’n  I  lu,  m 

dp +  It 


Equation  (2Ui  taken  at  the  boundaries  of  the  plasma  tat 
the  emitter  and  collector  sheath  interfaces)  forms  the 
plasma  boundary  conditions 
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is  written  as 
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where 

*,  .  ,>„/<»,«!  ti,M\ 

where  A  <  r )  is  the  ionization  coefficient  and  is  found  from 
consideration  of  ionization  kinetics  of  the  cesium  accord¬ 
ing  to  Lawless  [8|.  Its  solution  for  n  is 

«({)  -  B  sin  +  C)  (37) 

where  B  and  C arc  constants  of  integration  and  A  *  A  (  r  ). 
The  quantities  0„  and  0 , .  which  are  the  boundary  condi¬ 
tions  for  (37).  can  be  written  as  functions  of  r.  xt. 
and  Ay.: 

0»  »  fti (r.  \f.  x<-  d\J 

0,  »  0,(r.  xp.  Xi  •  ^X.)  (38) 

When  there  is  no  reflection.  0„and  0,  are  Kith  large,  i.e. 


°(e 


Significant  reflection  on  the  emitter  side  reduces  0,,  and  it 
may  indeed  attain  negative  values  for  sufficiently  strong 
reflection. 

The  density  equation  (35)  with  the  boundary  conditions 
0t)and  0,  is  a  linear  eigenvalue  problem:  its  solution  yields 
A  and  C  as  functions  of  0„  and  0,.  The  calculated  results 
are  shown  in  Fig.  3.  Since  A(  r )  is  a  function  of  r  from 
the  ionization  kinetics,  the  value  of  r  is  thus  determined 
by  a  function  of  0„  and  0,.  The  plasma  resistance  R  also 
can  be  expressed  in  terms  of  functions  of  0„  and  0,  through 
A  and  C  using  (27). 


I  (n,m  I  „ 


4  4/  sin  C 

__  —  j - |n 

s/5ff  A,  A 


•Pr) 

“(3) 


The  sheath  results  which  provide  j.  Q.  0„.  and  0,.  com¬ 
plete  the  isothermal  formulation.  The  sheath  theory  used 
is  an  exact  solution  to  the  Poisson  equation  and  collision- 
less  Boltzmann  equation  for  warm  plasma  ion  distribution 
with  a  10  percent  of  Bohm  speed  cutoff  velocity  to  ap¬ 
proximate  the  effect  of  the  collisional  presheath  |6|.  The 
results  are  summarized  below.  The  quantities  0„.  0,.  Q. 
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and  j  are  found  from  the  sheath  calculations  as  functions 
of  r.  Xa.  X< .  and  Ax,,  i.e.: 


ft.  *  ft.( T.  Xa-  Xc ■  Ax.) 

0.  «  0,(r.  xa-  Xc-  Ax.) 

Q  *  Q(r.  xa-  Xc-  Ax.) 
j  •  j(r.  Xa-  Xc-  Ax.). 

From  the  eigenvalue  problem  for  the  plasma  density  we 
then  hnd 


A(r)  -  0t).  (40) 


From  the  continuity  equation  for  current  we  find 


.  fsin(A  +  C)\ 

+  r  In  —  +  r  In  (  -  -  I 
On  W 

and  from  the  electron  momentum  equation  we  hnd 


(4! 


Xc  ”  Xa 


•  r  tn 


/sin  M  C )  \ 
\  sm  ( C )  / 


+  jR  + 


(42) 


These  three  previous  equations  determine  \t.  x, .  and  r 
when  A>.  is  given.  This  set  of  equations  is  valiJ  lor  jll 
A\,  Even  in  the  case  ot  A\,  S  0  when  there  is  no  reflec¬ 
tion.  the  calculations  differ  from  previous  isothermal  cal¬ 
culations  because  the  Boltzmann  assumption  on  the  elec¬ 
trons  ls  not  used  as  indicated  by  the  presence  of  a„  and 
<*.- 


Ul.  Calculated  Results  for  Ion  Reflection  and 
Trapped  Ions 

In  this  section  we  develop  isothermal  solutions  for  the 
thermionic  convener  with  the  emitter  sheath  phenomena 


.n.i 


TABLE  I 

Isothomai  Soli  tios  CosmTiov* 


CASE  1 

CASE  2 

Tc  -  IT  so  K 

Tc  «  l rso  K 

Tc  *  750  K 

T.~  -  7J0  K 

pf.  =  1  ton 

p„  »  1  torr 

d  *  10  mil 

d  *  10  nui 

4t  «  2.12  «V 

AC  -  2.07  eV 

4,:  m  1  GO  «V 

4r  »  1  60  .V 

Jn  =  20  imp 'em’ 

Jn  *  7.57anip/cm' 

J..+  •  1  JOxlO'4  imp/cm’ 

»  2. 10*10” 3  amp  cm5 

of  ion  reflection,  trapped  ions,  and  surface  emission  ions 
included.  Emitter  sheath  effects  on  thermionic  convener 
performance  can  be  divided  into  twe  categories:  I) 
changes  in  net  ion  flux  rate  into  the  sheath  which  affect 
plasma  density  directly  ,  and  2)  changes  in  sheath  poten¬ 
tial  distribution  which  affect  the  exchange  of  "hot" 
plasma  electrons  for  "cold"  emitter  ions  directly.  A  de¬ 
creased  influx  of  ions  into  the  sheath,  which  occurs  for  all 
three  emitter  sheath  phenomena,  increases  the  plasma 
density  at  the  neutral  plasma  emitter  sheath  interface. 
Theoretical  intuition  suggests  that  an  increased  plasma 
density  at  the  emitter  would  benefit  performance  by  re¬ 
ducing  resistance  through  the  plasma  and  therefore  reduc¬ 
ing  arc-drop.  However,  this  is  not  the  case.  While  the 
plasma  density  at  the  emitter  increases  slightly,  plasma 
density  at  the  collector  decreases.  Consequently,  total  re¬ 
sistance  increases. 

All  three  of  these  phenomena  increase  in  significance 
as  net  currem  density  through  the  convener  is  reduced. 
Each  of  these  reduces  the  net  ion  loss  rate  to  the  emitter 
and  consequently  increases  arc-drop  (therefore,  degrading 
performance  at  low  current  densities).  This  increase  in 
arc-drop  is  in  agreement  with  the  same  tendency  in  the 
experimental  tesults.  However,  the  experimental  results 
also  show  a  plateau  (of  low  arc -drop)  at  low  current  den¬ 
sity.  This  plateau  occurs  at  a  current  density  correspond¬ 
ing  to  significant  surface  ion  emission  and  is  therefore 
thought  to  occur  as  surface  emission  replaces  volume  ion- 
i/.mon  av  the  dominant  source  of  plasma  ions.  Unfortu¬ 
nately.  the  theoretical  calculations  cannot  be  carried  into 
this  region  because  the  collisionless  collector  sheai.i 
matching  (to  the  neutral  plasma)  fails. 

To  provide  a  realistic  framework  for  presenting  the  re¬ 
sults.  we  consider  the  convener  conditions  shown  as  case 

1  in  Table  I.  Case  2  is  shown  because  it  has  the  largest 
surface  emission  of  any  typical  thermionic  converter  op¬ 
erating  condition  (because  the  work  function  is  high  and 
the  temperature  is  also  high).  Instead  of  presenting  case 

2  separately,  we  demonstrate  the  effects  of  surface  emis- 
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sion  in  case  1  '  y  increasing  the  surface  emission  by  a 
factor  of  100  thereby  bringing  it  up  to  the  level  in  case  2. 
The  net  current  density  at  which  surface  emission  be¬ 
comes  significant  can  be  estimated  by  multiplying  by 
the  square  root  of  the  ion  to  electron  mass  ratio  (approx¬ 
imately  500).  In  case  1.  this  means  that  surface  emission 
becomes  significant  at  J  -  0.01  A/cm:  while  in  case  2 
significant  surface  emission  begins  at  J  -  1.0  A/cm:. 

IV.  Effects  of  Ion  Rfflection 

In  this  section  we  discuss  the  isothermal  results  for  case 
I  with  ion  reflection,  but  without  trapped  ions  and  with 
the  small  amount  of  surface  emission  ions  of  case  1 .  Fig. 
4  is  the  CV  diagram  for  this  case. 

The  dotted  line  extending  upward  from  point  A  is  the 
single  electron-. spelling  emitter  sheath  solution  How¬ 
ever.  we  have  not  taken  recombination  or  the  Schottkv 
effect  into  account  in  this  isothermal  formulation  which 
are  expected  to  become  important  at  current  densities  near 
J„.  The  interest  of  this  paper  begins  at  point  A.  where  the 
single  sheath  doubles  over.  Between  points  A  and  B. 
where  the  back  sheath  height  Ax  is  less  than  the  sheath 
height  Xf.  the  emitter  sheath  is  nonreflecting.  In  this  re¬ 
gion  the  sheath  heights  xr  and  Xc  remain  constant  while 
the  plasma  density  is  proportional  to  net  current  J  (the 
normalized  plasma  density  nc/J  is  constant).  Only  the 
back  sheath  height  Ax  changes  and  the  CV  curve  in  this 
region  is  Boltzmann  (the  arc-drop  is  constant).  Beginning 
at  point  B  and  continuing  to  point  C.  the  double  emitter 
sheath  reflects  plasma  ions  because  the  back  sheath  is 
larger  than  the  front  sheath:  in  other  words,  the  reflective 
potential  Ax,  m  Ax  ~  Xe  is  positive.  The  result  is  that 
net  ion  loss  rate  into  the  sheath  u  decreases  and  that  arc- 
drop  increases.  The  quantity  5  is  defined  as  the  mean  ion 
velocity  into  the  sheath  normalized  by  the  Bohm  speed 
yfkfjM.  The  dotted  curve  BD  is  the  same  double  sheath 
except  that  it  assumes  no  ions  are  reflected;  therefore,  5  is 
constant  and  arc -drop  is  constant.  The  two  curves  BC  and 
BD  are  almost  indistinguishable  because  the  increase  in 
arc -drop  is  small  until  the  net  current  density  is  extremely 
small.  The  reason  for  this  is  that  the  shift  speed  is  ap¬ 
proximately  u,  *  2.  and.  therefore,  a  large  increase  in 
reflective  potential  is  required  to  change  u  significantly 
(the  half- reflection  point  is  Ax,  *  4.0  or  approximately 
J  •  J*  exp  ( -4)  •  0.4  A/cm:).  The  shift  speed  u,  is 
defined  as  the  velocity  at  the  peak  of  the  incoming  ion 
distribution  again  normalized  by  the  Bohm  speed. 

The  curve  EF  is  the  single  electron-repelling  emitter 
sheath  case.  It  is  the  limiting  case  for  large  amounts  of 
trapped  ions  in  which  the  double  sheath  peak  has  been 
completely  suppressed  by  the  trapped  ions.  For  this  case, 
the  emitter  sheath  solutions  gives  u,  *  06.  This  curve  is 
not  topologically  connected  to  the  curve  ABC\  it  will  be 
shown  in  Section  V  that  trapped  ions  move  ABC  toward 
the  single  ion-repelling  sheath  case.  The  curve  is  much 
steeper  (a  faster  increase  in  arc-drop)  in  this  case  because 
u,  «  0  (the  half-point  in  ion  reflection  is  approximately  J 
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Fig.  5  Normalized  plasma  density  with  reflection 


*  8  A/cm:).  Curve  EG  is  the  single  ion-repelling  case 
assuming  no  reflection  and  is  therefore  a  Boltzmann  line 
with  constant  arc-drop. 

At  points  F"and  C  the  solutions  fail  at  the  collector.  The 
explanation  for  this  failure  is  best  given  by  examining 
Figs.  5-8. 

Fig.  5  is  the  normalized  plasma  density  through  the 
convener  gap.  The  highest  curve  with  no  reflection  Ax, 
«  0  has  the  largest  plasma  density  at  the  collector  but  the 
lowest  plasma  density  at  the  emitter.  Ion  reflection,  which 
decreases  the  ion  loss  rate  to  the  emitter,  raises  the  plasma 
density  at  the  emitter  but  lowers  the  plasma  density  at  the 
collector.  The  lower  plasma  density  at  the  collector  forces 
a  smaller  collector  sheath  heignt  to  pass  the  net  current 
density.  This  can  be  seen  from  (10).  Fig.  6  is  the  potential 
through  the  convener  under  the  same  reflection  conditions 
as  in  Fig.  5.  In  Fig.  6  the  first  two  spaces  on  the  left  make 
up  the  double  emitter  sheath,  and  the  last  space  on  the 
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Fig.  6  Potcantl  distribution  in  the  convener 
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Ft*  7.  Collector  sheath  failure. 

right  is  the  collector  sheath.  The  region  between  the  two 
sbeaths  is  the  neutral  plasma  region.  In  the  no-reflection 
case,  it  can  be  seen  that  the  potential  has  a  pronounced 
well  in  the  middle  This  is  the  result  of  the  large  plasma 
density  in  the  middle.  As  reflection  increases,  this  well 
disappears  on  the  collector  side  of  the  plasma  because  re¬ 
sistive  drop  f*"tre  (due  low  plasma  density)  increases 
to  the  degree  -%t  it  is  greater  than  the  ambipolar  rise  (due 
to  decre*  i  »  iensity  toward  the  collector).  Simulta¬ 
neously  *v  .....  fs.^sma  potential  gradient  at  the  collector  be¬ 
coming  i.ega';’v«,  ie  collector  sheath  goes  toward  zero 
height.  Fig.  "  ^.-ws  the  critical  collector  sheath  quan¬ 
tities  as  the  collector  sheath  failure  occurs.  Collector 
sheath  height  xc  goes  toward  zero,  the  shift  speed  u,e  goes 
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Fig.  I.  Emitter  theath  during  reflection. 

toward  negative  infinity,  and  the  ion  loss  rate  to  the  col¬ 
lector  ue  is  driven  to  zero.  The  two  preceding  quantities 
ue  and  u,e  are  defined  at  the  collector  sheath  as  u  and  /<, 
were  at  the  emitter  sheath.  Fig.  8  shows  the  changes  in 
the  emitter  sheath  height,  ion  shift  speed  and  ion  loss  rate. 
When  the  collector  sheath  failure  occurs,  the  ion  loss  rate 
to  the  collector  is  zero  ( 3,  *  0)  and  the  corresponding 
plasma  ion  distribution  at  the  collector  is  bunched  at  zero 
velocity  (u,e  «  -on).  While  the  mathematics  hold  self- 
consistantly  until  5r  «  0,  the  physics  is  clearly  poor  at 
this  point  because  5r  ■  0  demands  that  the  plasma  ions  at 
the  collector  have  zero  energy  (zero  temperature  and  zero 
mean  velocity).  An  estimate  of  when  the  physics  becomes 
poor  is  u,r  «  0.  At  this  point  the  net  ion  loss  rate  is  close 
to  the  thermal  speed.  A  second  physical  difficulty  that  oc¬ 
curs  with  collector  sheath  failure  is  that  the  electron  Mach 
number  there  Qc  (from  (10))  becomes 


because  the  collector  sheath  height  approaches  zero  ( ac¬ 
tually  about  0.001 ).  In  the  present  continuum  formulation 
of  the  plasma  region,  it  was  assumed  in  (13)  that  Q ,  is 
small  so  that  the  electron  momentum  term  urdu,/dx  can 
be  neglected. 

One  could  take  the  solution  below  the  collector  sheath 
failure  point  if  2r  could  attain  negative  values  or  if  Q , 
could  attain  values  larger  than  ■JT/ir.  There  is  no  physical 
basis  for  assuming  that  ur  can  become  negative  since  the 
collector  emits  nothing.  However,  there  is  a  physical  ba¬ 
sis  for  allowing  Q \  to  be  larger  than  JTfx  (an  electron 
distribution  shift)  as  can  be  seen  in  Fig.  6:  the  potential 
drop  nearing  the  collector  becomes  progressively  more 
electron  accelerating  as  the  collector  sheath  fails,  and. 
therefore,  the  electron  distribution  should  be  shifted  as 
the  ion  distribution  is  in  an  electron-repelling  sheath. 
However,  this  would  clearly  invalidate  the  assumption 
that  the  electron  momentum  term  is  negligible.  Therefore, 
the  momentum  term  must  be  added  to  explore  further  in 
this  direction  and  this  has  not  been  done  because  of  the 
resulting  complexity  in  the  equations. 
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Comparison  of  Fig.  7  to  Fig.  8  at  the  collector  sheath 
failure  point  (Ax,  *  2.5,  uc  *  0)  shows  that  the  ion  loss 
rate  to  the  emitter  is  positive.  At  this  point  the  plasma  is 
still  ignited  and  generating  ions  as  can  be  seen  from  Figs. 
9  and  10.  The  ionization  coefficient  A  has  dropped  by  50 
percent,  but  the  plasma  electron  temperature  has  dropped 
by  only  5  percent.  Finally,  we  note  in  Fig.  11  that  the 
normalized  plasma  resistance  R  has  risen  by  almost  100 
percent.  This  is  responsible  for  the  increase  in  arc-drop 
and  the  decrease  in  performance.  Plasma  resistance  in¬ 
creases  in  response  to  reflection  because  the  loss  of  plasma 
electron  energy  to  the  emitter  is  more  important  than  the 
loss  of  ionization  energy  to  the  emitter.  Ion  reflection  at 
the  emitter  increases  the  normalized  plasma  density  there, 
and  consequently  increases  the  normalized  loss  of  plasma 
electron  energy  there.  The  basis  of  this  can  be  seen  from 
conservation  of  electron  energy  (22): 

r  *  I  ”  \)V<  -  \j.Vf.  (43) 

The  ion  energy  loss  term  is  generally  small  compared  to 
the  electron  energy  loss  term: 


Fig.  10.  Plum*  election  temperature. 


not  see  any  plateau  or  decrease  in  arc-drop  as  net  -urrent 
density  is  decreased  in  the  present  calculations. 


\i<  vt< 

Wj  ~~J  v“ 


0(0.02). 


(44) 


Therefore,  we  take  the  electron  energy  equation  as 


r-I-i  jV4.  (45) 

Since  r  is  nearly  constant  (because  of  the  ionization  ki¬ 
netics),  the  product/)^  is  nearly  constant.  Ion  reflection 
decreases  j  (because  the  normalized  plasma  density  in¬ 
creases)  and  therefore  increases  arc-drop  Vt  ( makes  V4  a 
more  negative  number). 

If  the  equations  are  reformulated  in  such  a  way  as  to  be 
valid  past  the  collector  sheath  failure  point,  then  we  can 
eventually  expect  to  see  a  decrease  in  arc-drop  and  a  low- 
current  plateau  as  the  electron  temperature  approaches  1 
(the  ignited  plasma  is  extinguished  and  the  ionization 
source  is  surface  emission).  This  can  be  seen  from  (43). 
However,  as  we  see,  the  collector  failure  occurs  before  r 
has  dropped  more  than  5  percent.  Consequently,  we  do 


V.  Effects  of  Trapped  Ions 

Fig.  12  shows  the  effect  of  trapped  ions  on  the  CFi 1  ar- 
acteristics.  In  this  section  the  trapped  ion  distribution  is 
assumed  to  have  the  temperature  of  plasma  ion  distribu¬ 
tion,  and  100-percent  trapped  ions  (  fir  ■  1.0)  is  defined 
to  complete  the  ion  distribution  at  the  double  emitter 
sheath  peak  such  that  one  has  a  Maxwellian  distribution 
there.  Based  on  physical  reasoning  about  the  trapping 
mechanism,  one  expects  on  the  order  of  10.  Also,  some 
trapping  calculations  have  been  done  for  approximate 
sheath  formulations  [9],  [10]  which  support  this. 

Curve  AHIJ  is  the  CV  characteristic  for  f„  ■  0.10.  At 
point  A  there  cannot  be  any  trapped  ions  since  the  back 
sheath  height  Ax  is  zero.  Therefore,  the  trapped  CV 
merges  into  the  nontrepped  curve  there.  The  actual  amount 
of  trapped  ions  on  the/,r  **  0. 10  curve  increases  from  zero 
at  point  A  to  the  full  10  percent  of  a  thermal  distribution 
at  point  H  where  the  back  sheath  height  Ax  is  equal  to  the 
sheath  height  xe-  The  shift  speed  increases  on  AH  from 
1.95  to  3.00.  This  corresponds  to  what  is  seen  in  Fig.  12 
where  Ax  <  Xe- 
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The  rise  in  shift  speed  has  been  limited  to  3.00.  This 
limit  »  placed  on  the  shift  speed  because  a  sheath  with 
height  of  about  1.0  should  not  have  a  presheath  region 
capable  of  shifting  the  entire  distribution  so  far.  In  fact, 
limiting  the  shift  speed  is  equivalent  to  increasing  the  cut¬ 
off  speed  for  the  ion  distribution  function. 

The  arc-diop  decreases  as  a  result  of  the  increase  in  u, 
and  the  consequent  increase  in  the  net  ion  loss  rate  to  the 
emitter.  A  “hump"  can  be  seen  on  AH  where  the  shift 
speed  hits  3.00.  The  arc-drop  is  lowest  on  this  "hump" 
because  the  shift  speed  is  at  its  maximum  of  3.00.  Be¬ 
tween  points  H  and  /  the  back  sheath  height  remains  equal 
to  the  sheath  height.  Ax  -  Xe  *  X.  *  0.  On  this  segment. 
«,  decreases  to  1.25.  therefore  increasing  arc-drop. 

From  point  /  to  point  J,  the  shift  speed  remains  constant 
at  1 .25  and  the  ion  loss  rate  decreases  because  of  reflec¬ 
tion.  The  other  trapped  cases  /„  *  0.2.  0.3.  and  0.4  have 
not  been  connected  because  they  hit  the  3.00  maximum 
shift  speed  much  sooner  than  in  the /,  *  0. 1  case. 

Point  J  is  the  collector  sheath  failure  point.  Each  of  the 
f„  »  0.2.  0.3.  and  0.4  curves  begins  at  Ax,  *  0  and  ends 
at  the  collector  sheath  failure  point.  It  should  be  noted  that 
each  of  the  trapped  ion  curves  fails  at  a  higher  current 
than  the  last  because  the  shift  speed  is  lower. 

VI.  Effects  of  Emoter  Surface  Emission 

Fig.  12  shows  the  effect  of  surface  emission  on  the  /, 
«  0.10  curve:  surface  emission  is  added  by  multiplying 
the  actual  small  amount  of  surface  emission  in  case  1  by 
a  factor  of  100.  This  brings  the  surface  emission  up  to  the 
level  in  case  2.  malting  it  significant  at  J  •  1.0  A/cm:. 
It  can  hr.  seen  that  surface  emission  increases  arc-drop:  it 
does  so  in  exactly  the  same  way  as  reflection  or  trapped 
ions  do— it  decreases  the  net  loss  rate  of  ions  to  the 
emitter. 

VQ.  Comparison  wtth  Experimental  Results  and 
Conclusions 

Fig.  13  superimposes  the  isothermal  results  of  Fig.  12 
on  the  experimental  results  for  a  cesium  reservoir  tem¬ 
perature  of  SSI  1C  which  produces  a  l-torr  neutral  cesium 


Fig.  13.  Isothermal  anti  experimental  Ct  diagrams 

pressure.  The  experimental  results  are  from  |11].  The 
point  of  this  comparison  is  that  the  steepness  of  the  CV 
characteristic  in  the  experimental  converters  can  be  ex¬ 
plained  by  a  decreasing  ion  loss  rate  to  the  emitter.  We 
have  shown  that  all  three  of  the  expected  emitter  sheath 
phenomena  decrease  the  ion  loss  rate  to  the  emitter.  We 
cannot  calculate  the  amount  of  trapped  ions  in  a  collision¬ 
less  sheath  without  knowledge  of  the  coliisional  pro¬ 
cesses.  However,  the  experimental  CV  suggests  that  if  the 
amount  of  trapped  ions  ( /,, )  increases  from  0  percent  at 
J  «  14  A/cm:  (the  double  sheath  formation  point)  to  10 
percent  at  J  •  2  A/cn r.  then  the  steepness  could  result 
from  trapped  ions  reducing  the  ion  loss  rate  to  the  emitter. 
Since  these  percentages  are  based  on  a  thermal  distribu¬ 
tion- of  ions,  they  seem  physically  reasonable.  Unfortu¬ 
nately.  the  collector  sheath  failure  prevents  us  from  going 
to  the  point  in  the  calculations  where  r  drops  enough  to 
make  surface  emission  the  source  of  ions. 

The  experimental  curve  is  nearly  a  constant  0.05  V  be¬ 
low  the  isothermal  result  ( f„  »  0. 10)  except  at  high  cur¬ 
rent  densities  and  at  the  "hump."  Comparison  of  the 
curves  at  high  current  density  is  not  valid  since  neither 
the  Schottky  effect  nor  recombination  has  been  included. 
The  Schottky  effect  is  important  above  12  A/cm:  in  this 
case  because  the  emitter  sheath  is  single  electron  repel¬ 
ling  (to  the  plasma)  and  therefore  puts  a  strong  electric 
held  against  the  emitter  with  the  appropriate  sign.  Recom¬ 
bination  is  also  potentially  important  because  the  plasma 
density  scales  with  current  density,  and  at  high  current 
densities  the  plasma  density  in  the  middle  of  the  convener 
approaches  the  Saha  density.  The  0.05-V  difference  may 
or  may  not  be  explained  by  a  discrepency  in  the  assumed 
collector  work  function.  At  750  K  the  collector  emits  es¬ 
sentially  nothing  and  therefore  any  change  in  the  collector 
work  function  directly  affects  output  voltage.  If  the  col¬ 
lector  work  function  were  in  fact  1.65  instead  of  1 .60  V, 
then  the  isothermal  result  would  lie  nearly  on  top  of  the 
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experimental  result.  We  have  not  adjusted  the  assumed 
collector  work  function  so  as  to  illustrate  the  importance 
of  it  and  therefore  the  importance  of  the  surface  physics 
of  the  adsorbed  cesium  layer.  The  “hump”  should  not  be 
taken  as  an  expected  experimental  result  since  it  results 
from  the  interaction  of  the  trapped  ions  with  the  plasma- 
emitter  sheath  interface.  Instead  it  should  be  taken  as  a 
second  reason  (in  addition  to  the  cutoff  of  the  ion  distri¬ 
bution)  for  further  study  of  the  matching  region  between 
the  collisionless  sheath  and  the  neutral  plasma. 
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ASYMPTOTIC  SHEATH  THEORY 


The  Asymptotic  Sheath  Theory  developed  under  this  grant  is  set  forth  in  the  following 
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Few  exact  solutions  for  collisional  presheaths  exist  because  of  the  difficulty  of  simultaneously 
satisfying  both  the  collisional  Boltzmann  equation  and  the  Poisson  equation.  The  exact 
solutions  that  do  exist  are  for  very  specialized  collision  terms  such  as  constant  cross-section 
charge  exchange  with  cold  neutrals.  The  present  paper  presents  an  asymptotic  method  which  is 
applicable  to  a  variety  of  collision  terms  and  is  applied  in  particular  to  constant  collision 
frequency  charge  exchange  with  noncold  neutrals.  Constant  collision  frequency  and  constant 
cross-section  collision  with  cold  neutral  results  are  also  presented.  The  first-order  terms  for  the 
presheath  potential  rise  and  ion  distribution  functions  are  calculated  and  it  is  shown  that 
second-  and  higher-order  terms  can  be  calculated  using  a  multiexponential  expansion  for 
presheath  potential  rise.  The  first-order  cold  neutral  constant  cross-section  results  correspond 
well  to  the  exact  solution.  The  calculated  presheath  potential  rises  are  of  the  order  expected 
from  the  Bohm  criterion,  and  in  some  of  the  specialized  cold  neutral  cases,  exactly  kT,/2.  The 
presheath  potential  rise  is  reduced  by  a  neutral  plasma  potential  gradient  which  accelerates 
ions  toward  the  presheath.  In  all  cases  the  collisional  presheath  is  asymptotically  matched  to 
both  the  neutral  plasma  and  the  collisionless  sheath. 


L  INTRODUCTION 

The  majority  of  plasma-surface  interaction  work 
matches  a  neutral  plasma  to  a  collisionless  sheath  without 
detailed  consideration  of  a  collisional  presheath.  However, 
the  collisional  presheath  structure  is  of  great  interest.  Sheath 
theory,  beginning  with  Bohm, 1  tends  to  assume  that  the  plas¬ 
ma  ion  distribution  is  cold  so  that  a  minimum  presheath 
potential  rise  may  be  calculated,  which  makes  the  collision¬ 
less  sheath  self-consistent  Harrison  and  Thompson2  genera¬ 
lize  the  Bohm  criterion  to  noncold  km  distributions;  how¬ 
ever,  the  result  is  sensitive  to  the  density  of  the  low  energy 
tail  of  the  km  distribution,  which  in  turn  is  strongly  affected 
by  the  collisional  presheath.  And,  a  second  difficulty  in  the 
absence  of  a  collisional  presheath  is  that  the  collisionless 
sheath  and  the  surface  beyond  it  may  return  no  ions  or  a 
''.on  thermal  distribution  of  ions  which  the  collisional  pre- 
sheath  must  match  to  the  neutral  plasma  region. 

Some  exact  solutions  exist  for  presheaths;  notable  is  the 
work  of  Ecker  and  Kanne3  and  Riemann,4  who  derive  exact 
solutions  for  collision  terms  based  on  charge  exchange  with 
cold  neutrals  and  Emmert  et  a/.,3  who  derive  an  exact  colli¬ 
sionless  solution  in  which  there  is  an  ionization  source.  In 
the  present  paper  an  asymptotically  correct  collisional  pre¬ 
sheath  theory  is  developed  which  can  be  applied  to  a  less 
restrictive  range  of  collision  terms.  Potential  in  the  pre¬ 
sheath  is  expanded  as  a  multiexponential  series  and  the  dis¬ 
tribution  functions  are  expanded  in  terms  of  presheath  po¬ 
tential  rise.  First-order  approximations  are  calculated  for 
both  constant  collision  frequency  and  constant  cross-section 
charge  exchange  collisions. 

IL  FIRST-ORDER  ASYMPTOTIC  POTENTIAL 
FORMULATION 

In  this  section  it  is  assumed  that  the  potential  in  the 
collisional  presheath  is  of  the  form 

Um  U0  +  At/  —  ax  +  *A,f  (1) 


where  t/0  *  ax  is  the  assumed  linear  potential  in  the  neutral 
plasma  and  A  (7  ■»  is  the  additional  potential  rise  in  the 
collisional  presheath,  as  shown  in  Fig.  1.  In  this  paper  the 
convention  used  is  that  U  *  qi,  where  q  is  the  electron 
charge  and  l  is  potential  in  electron  volts  so  that  U  has  units 
of  energy.  In  addition,  potential  is  defined  in  the  reverse  of 
the  usual  sign  convention  so  that  increasing  potential  repels 
electrons.  With  these  conventions,  the  Boltzmann  equation 
can  be  written  as 


In  Eq.  (2)  and  those  following,  the  ±  denotes  the  sign  of 
the  charged  species  in  question;  the  upper  sign  refers  to  posi¬ 
tively  charged  ions  and  the  lower  sign  to  electrons.  The 
Boltzmann  equation  is  expressed  in  terms  of  A  U,  which  will 
be  the  expansion  variable  in  the  presheath: 

V0*U'£ r/±~(/*A£/+a)lr,B(:r)  •  (3) 

<?AI/  m  dv  \dt  )t 

The  distribution  function  is  then  expanded  as 

/-/o(v)  +  ACft(t>)+AJ/j(i>)-l---.  (4) 

so  that  the  derivatives  are 

^-/i(»)  +  2A£//2(u)  +  3At/,/,(u)+  ••  (5) 

and 

|f«%(tf)  +  AC/^(u)  +  Al/J^-(u)  +  ••  (6) 

do  dv  do  dv 

Substitution  of  ( 5 )  and  (6)  into  the  Boltzmann  equation  ( 3 ) 
yields  the  terms 
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(7a) 


I: 


AUi 


vPUv)  ±S.§k(V)±SL^ I („) « 

m  dv  m  dv 


AU3z  2c@f2(.v)  ±  —  (p)  ±  —  (t>) 

m  av  m  av 


(7b) 

(7c) 


AIT:  mfif.(v)  ± 


£#Lj-(„)±.2.^L(„) 

m  av  m  av 


The  quantity  #  representing  the  presheath  potential  rise,  is  determined  from  the  Poisson  equation 


(7d) 


=  f*(J  MvAU)d»  —  f  /(c,A£/)di>j,  (8) 

where  q  is  the  electron  charge.  It  is  assumed  that  the  ions  are  singly  ionized  for  simplicity.  The  Poisson  equation  (8)  is 
expanded  as 

0*AU-  W[(j~  /,(»)</«>- J*  f.,{v)d^jAU  +  (J“  fa  (v)dv  —  J*"  /«(tf)4fo)A£f*+ •••],  (9) 

where  charge  neutrality  at  AU  **  0  has  eliminated  the  terms  containing /„  and /o : 


«o« 


{o)dv. 


The  quantity  n,  is  the  neutral  plasma  density  of  the  asymptotic  presheath,  not  of  the  neutral  plasma. 


HI.  FIRST-ORDER  SOLUTION  WITH  A  CONSTANT  COLLISION  FREQUENCY  CHARGE  EXCHANGE  COLLISION  TERM 
The  constant  collision  frequency  charge  exchange  collision  term  is  modeled  as 

(x)  ~'rn f(u)du-f(v)J  /„{u)duj,  (10) 

where/,  (o)  is  the  neutral  distribution  and  r  is  the  collision  time.  Previous  work  has  assumed  cold  neutrals  and  results  in  an  in¬ 
tegral  equation  which  is  solvable  only  for  constant  collision  cross  section.4 


A.  Zero  plasma  potential  gradient  (a*0) 

In  this  case  Eqs-  (7)  become 

1:  °*mK^,<0)  I  -/•<»>  J  /.(«)<*«). 

A  U:  vfffai")  +  —  *~(»)  (/„(«)  f  /n  (u)du -/t  (v)  f  fm(u)du), 

m  dv  mn\  J-.  J-m  J  (11) 

ACf':  «C(»)4%li(»)-U/,(r)  H  f»(u)du-fm{v)  f"  f„(u)du\ 
m  dv  m,  \  J  -  m  J- m  / 

Under  the  assumption  that  the  neutral  distribution  is  Maxwellian/,  (v)  -  n^m/lnkT  exp(  —  mif/lkT),  the  solution  to 
(11)  is 

fn (®)  *  (l/k7T/(((p),  fn> 


/■O')  -  ( UnkT)flim _„(»)- 

Thus 

/(rAtO-Ct,M"‘n/,(i»),  ^  (13) 

which  is  the  expected  result.  In  this  case  the  mean  ion  velocity  is  zero  throughout  the  collisions!  presheath  since  charge 
exchange  collisions  conserve  ions  and  the  mean  ion  velocity  in  the  neutral  plasma  is  zero.  Thus,  if  a  »  0,  constant  collision 
frequency  charge  exchange  collisions  do  not  shift  the  ion  distribution  upward  in  velocity.  This  presheath  can  be  matched  to  a 
collisionless  sheath  only  if  the  collisionless  sheath  returns  all  the  ions  entering  it  from  the  c  j’lisi ^nal  presheath. 

With  electron  density  assumed  to  follow 
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<KTr\  _  J-^V/kT,) 

ji,(AC/)  =  nQe  , 

the  Poisson  equation  (9)  yields,  to  first  order, 

01  =  4ir<i3no(l/kT+  1/kT.), 

which  is  the  length  scale  of  the  Debye  length.  Thus  for  a  =  0  the  collisional  presheath  is  not  distinct  from  the  collisionless 
sheath  since  there  is  no  separate  collisional  presheath  length  scale. 

B.  Nonzero  plasma  potential  gradient  (a^O) 

Under  this  condition  there  is  a  net  flux  of  ions  from  the  plasma  into  the  sheath,  which  allows  the  construction  of  a 
collisional  presheath  that  accelerates  the  ions  and  depopulates  the  ion  distribution  of  returning  ions.  Thus  the  collisional 
presheath  may  be  correctly  matched  to  the  collisionless  sheath  which  returns  no  ions.  In  this  case  ( 7a )  and  (7b)  can  be  written 

as 

2-&(v)~—[/Av)n0-nJ0(v)],  (14a) 

m  dv  m „ 

*/,(«;)  +  £&  +  -%-  00  =—[/.(»)«,  -«y,(»)].  04b) 

m  dv  m  dv  rnm 

The  solution  to  Eqs.  ( 14)  are 

, ,  ,  (  —mv\m  I  m  C  (  mu2  ,  mu\  , 

fcyliHT  J +„r“  <15) 

and 


ere 

»o*J  f0(v)dv 
i 

«i*J  ft(o)dv. 


GTopfsd)  _£&(«)!*,  + cl.  (16) 

ar/L  ar  \l  arrkT  \lkTJ  a  dv  I  J 


The  constant  of  integration  in  (15)  has  been  set  so  that  f0  goes  to  zero  at  —  ao;/0  goes  to  zero  at  oo  regardless  of  the 
constant  of  integration.  Equation  (17)  is  immediately  satisfied  by  (15).  The  constant  of  integration  C  in  (16)  must  beset  so 
that  (18),  which  represents  self-consistency,  is  satisfied.  It  can  be  seen  from  (16)  that/,  goes  to  zero  at  —  oo  and  oo  regardless 
of  the  constant  C.  From  ( 18 ),  then 


c  -  4  -  JT.  -K  -  -  5)  JT  + <  -  S?)  *  *] 

*  [  VI  “Kw)] ' ' + [/:.  -k  -  ^  if  ^ + =)*«■>*  *] 

The  exponential  sheath  rise)?  is  determined  from  the  Poisson  equation  under  the  simplifying  assumption  that 


One  might  expect  that  the  approximation  should  be  */t0exp(  —  U/kT,);  however,  this  cannot  be  true  in  the  asymptotic 
presheath  because  n,  must  approach  n0  as  U approaches  negative  infinity.  With  (20)  the  Poisson  equation  (9)  to  first  order 
becomes 

B1  *  4rgJ(n,  +  rif/kT, ).  (21) 

Since  the  ion  density  is  only  calculated  to  first  order,  the  same  will  be  done  for  the  electron  density  in  (20). 

To  obtain  a  particular  solution  it  is  assumed  here  that  the  collisionless  sheath  to  which  the  collisional  presheath  is  joined 
at  At/*  A U*  returns  no  ions.  In  particular, 

f  f[v,AV)dv-0,  (22) 
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J°  /0(u)<fu  + At/*  J°  /,(»)<&  *0 

t 

J"°  vf(v,&U*)dv  =  0, 

J°  vf0(v)dv  +  AC/*  J°  qf,(r)du*0. 


Berate  the  approximation  is  only  first  order,  it  is  not  possible  to  impose  the  condition  that  f(v)  is  uniformly  zero  for  returning 

Equations  ( 23 )  and  ( 25 )  represent  zero  returning  ion  density  and  zero  returning  ion  flux.  When  higher-order  terms  are 
iiacluded,  the  conditions  of  zero  returning  ion  momentum  flux,  zero  returning  ion  energy  flux,  etc.,  can  be  applied  in  succes¬ 
sor  Equations  ( 2 1 ) ,  ( 23 ) ,  and  ( 25 )  are  solved  for  n , ,  0,  and  At/*,  with  all  other  quantities  assumed  constant.  Equation  (21) 
immediately  satisfies  the  Bohm  criterion  at  At/ *  At/*  for  the  first-order  approximation 

*i  +  *o/kT,  >0.  (26) 

The  Poisson  equation  (21 )  can  be  written  as 

0UJO«1  +  kT.(nt/n0),  (27) 

where 

A.0  — >//c7’#/4irfin0  (28) 

is  the  Debye  length.  It  i:  expected  that  the  length  scale  of  the  presheath  should  be  of  the  order  0  —  1/A,,  where  A,  is  the  ion 
mean  free  path.  In  the  circumstance  that  the  Debye  length  is  small  compared  to  the  ion  mean  free  path,  the  product  0 2  A  is 
small  and 

a,  sx  —  rto/kT,.  (29) 

The  neutral  plasma  region  is  matched  to  the  collisional  presheath  also  at  At/  *  At/  *,  as  shown  in  Fig.  1,  to  produce  a 
three -scale  uniform  asymptotic  solution.  In  particular,  assuming  constant  collision  frequencies,  the  momentum  equations 
become 

(3o, 

\  rr  J  dx  dx  t 


(,„  "»,r? \dn  dU  m,  r, 

[ 

w,r  l\dn  dU  »"#r, 

\.kT--—rz- 


where 


.coHwcrW  worn 


/u«u,*au 


wra  piowna- 


F1G.  1.  Asymptotically  correct  potential  in  the  collisional  presheath. 
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(32) 


in  g  id!* 

T"*  ~  noP  ’ » 
dx  kT, 


—  =  a  +  0MJ*.  (34) 

dx 

jh«  quantity  n  is  the  plasma  density  at  the  matching  point  L.U*  and  I~,  and  I~,  are,  respectively,  the  ion  and  electron  net  fluxes. 
Nondimensionalization  results  in 

A  *  (ar/m)^m/lkT ,  (35) 

B=0kT/a,  (36) 

R,  =  T,/T,  (37) 

a  s , jm/lkT  v,  (38) 

whenAtndR,  are  the  parameters  and  2)  is  a  function  of  A  and  A,. The  quantity  .4  represents  the  nondimensional  asymptotic 
presheath  potential  gradient,  B  represents  the  nondimensional  exponential  presheath  rise,  R,  is  the  electron  to  neutral 
temperature  ratio,  and  a  is  the  nondimensional  velocity.  The  distribution  functions  can  then  be  written  as 

- r  exp(-{’  +  4)«  (39) 

n0yJm/2kT  -JvA 


Ft(a><44) 


(no/kT,  Um/2kT 


-  R.B  [-Up!  -  S ’>  -  ^  -  i)  £.  “■>(  -  V1  +  *>■])*  +  4 

_  _  ^  r .  "p(  ■  ^ )  r  ^  f  ^ 

+ ^  r.  “K  -  ^  -  t)  r  ■ + 

x  [i.  «p(  -  c ■>  -  ±  “K  -  4-)  £ .  «<  -  n‘  *  ' 


Thus  (23)  and  (25)  become 


(aAJi)do)  =■  0 


mF0(uA  )do) 


(aAJ)  da  “  0. 


Figure  2  presents  the  presheath  potential  rise  W/kT,  and  the  nondimensional  exponential  rise  B  as  a  function  of  the 
nondimensional  asymptotic  presheath  potential  gradient  A  for  a  range  of  electron  to  neutral  temperature  ratios  R, .  As  would 
be  intuitively  expected,  the  presheath  potential  rise  decreases  with  increasing  A.  Figure  3  presents  the  ion  distribution 
functions  at  the  neutral  plasma-collisional  presheath  interface  F0(a),  the  first-order  correction  to  the  distribution  function 
f i  («),  and  the  resulting  distribution  function  at  the  collisional  presheath-collisionless  sheath  interface  F0{a)  a-  &U*Fx(a). 
Although  the  resulting  distribution  is  not  uniformly  zero  for  a  <  0,  its  net  returning  density  and  flux  are  zero  by  (42 )  and 
( 43 ) .  It  is  expected  that  higher-order  corrections  to  the  distribution  function  and  potential  with  the  corresponding  application 
of  higher-order  moment  conditions  of  zero  returning  momentum,  energy,  etc.,  will  converge  the  returning  distribution 
function  toward  a  uniform  zero. 

In  the  limit  of  cold  neutrals,  the  constant  collision  frequency  charge  exchange  solution  is  considerably  simplified.  Equa¬ 
tions  (14a)  and  (14b)  become 
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—  —  (»)  *  — — [M.  (v)n0  -  n  J0{v)  ] 
m  dv  m. 


FIG  2-  Constant  collision  frequency 
presheath  potential  rise. 


tf/,(r>+  — +  — 4r(t,)  « <»)«,- «*/i(»)J- 

m  dv  m  dv  rnm 

s  solutions  to  (44)  and  (45)  are 

f#i0(m/ar)exp(  -  mu/ar),  o>0, 
fo\v>  *  |q  V<0, 


„pf.&2i-=;)(c-).  »<o, 

V  2a  ar/ 


•TOO  -525  -ISO  -l  75  000  >75  350  525 


5  fttys  Ru«.  VqI.3Q.no.  6.  Jaw  1M7 


FIG.  3.  Constant  collision  frequency  ion  distri¬ 
butions  in  the  neutral  plasma  and  at  the  pre- 
tbeatb-sbeath  boundary. 
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(48) 


such  that 

C*  —  C~  -  (m/ar)[«,  —  (fi/a)n0 ]. 

Equation  (46)  immediately  satisfies  n0  —  /"  m/0(v)dv.  No  returning  ions  implies  that 
C~  —  0 


C*  *  (m/ar)[nl  -  (0/a)rro] 

since y, on  r  < 0 is  already  zero.  The  final  condition  is  then  that  n,  =  f£f,(v)dv,  or 

-  f‘  op(  -  ^  ^  +  ^-2- Y  f  exp(  1*- 

Jo  V  2a  arJ  l  or  a  or  a  \arJ  Jo  \  2a  /  1 

The  application  of  it,  *  —  H^/kT,  yields 

4*4_/,  f_y  .uHr.,,/  ^  ,  ,if,  f*.., 

\  2m  3 


In  this  case,  A  U*  is  defined  by 
/0((r)  +  A(/*/,((T)-0. 

which  yidds 

W/kT.  -  \/(f}kT,/a  +  1), 
as  expected.  In  the  limit  of  pat3  /2m  — 0  we  have 
fikT./a  *  l 


(49) 

(50) 

(51) 

(52) 

(53) 

(54) 

(55) 


W-kT./Z  (56) 

which  corresponds  to  the  Bohm  criterion.  Figure  4  presents  the  variation  of  B  ~  0kT,/a,  with  0ar*/2m  for  the  cold  neutral 
case.  A  particular/?  for  the  parameters  can  be  conveniently  found  by  drawing  a  line  from  the  origin,  with  slope  Imk  T,/i3a1, 
so  that  the  intersection  is  the  solution.  Figure  S  presents  an  example  cold  neutral  ion  distribution.  Examination  of  the  ion 
distribution  function  at  o  ■>  0  shows  that  the  slope  is  discontinuous.  This  is  because  the  neutral  source  is  a  delta  function  at 
v  *  0.  It  appears  that  the  Bohm  criterion  cannot  be  satisfied  at  A  U  •  because  the  integral  Joi/X^/v2  ]dv  is  singular,  however, 
the  use  of  this  integral  in  the  Bohm  criterion  assumes  that  the  ions  accelerated  are  not  replaced.  In  this  case  the  ions 
accelerated  hum  u  m  0  are  replaced  by  ions  from  the  cold  neutral  distribution  which,  of  course,  is  a  delta  function  at  v  *  0. 


IV.  FIRST-ORDER  SOLUTION  WITH  A  OUASICONSTANT  CROSS-SECTION  COLLISION  TERM 
Fust-order  asymptotic  solutions  can  also  be  developed  for  a  quasiconstant  cross-section  collision  term 

|»-a»|dai-J  .  (57) 


FIG.  4.  Constant  coUistoo  frequency  pmheath 
nse  with  cold  neutrals. 
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FIG.  S.  Constant  collision  frequency  ion  distri¬ 
bution  with  cold  neutrals. 


This  collision  term  is  not  really  constant  cross  section  because  it  is  a  one-dimensional  representation  which  does  not  take  into 
account  average  velocities  in  the  other  two  dimensions.  However,  this  collision  term  corresponds  to  that  commonly  called 
constant  cross  section.  The  application  of  this  term  leads  to  a  set  of  integro-diflerential  equations  which  can  be  at  least 
approximately  solved,  and  in  the  cold  neutral  case  it  leads  to  readily  soluble  first-order  differential  equations.  The  cold  neutral 
case  presented  here  corresponds  to  that  which  can  be  solved  exactly  ( Riemann4 ) .  Unfortunately,  though,  the  exact  solution 
method  is  not  extensible  to  noncold  neutrals.  The  cold  neutral  collision  term  is 


(Jrl  ftu)\u\du-of{v)nm\v\ 

and  the  zero-order  Boltzmann  equation  term  ( 7a)  becomes 

»5(|>)J  f0(u)\u\du  —  an„  |u|/o(c). 


»  cm. 


m  dv 

for  which  the  solution  is 


/o(°) 


fi  jamn.  (  <mnn  \ 

"%/ta/—  "n-ir'7  c>0’ 


o. 


The  first-order  Boltzmann  term  is 

q6T/,(o)  +  ~~(c)  +  — ^L(tf 
m  dv  m  dv 

for  which  the  solution  is 

f«p[ 


r  <C. 


\u\du  -anK ]i>|/,(t>). 


/,(e) 


The  jump  oondition  itv«0  must  be  satisfied  in  ( 61 ) : 

No  returning  ions,  C'»0,  and  the  application  of  ( 63 )  to  ( 62 )  yields 
C *  —  no(0 /a)\2/Vyfcrmnm/a . 

The  oollisional  presheath-collisionless  sheath  boundary  &U*  is  again 


v  <0. 


(58) 

(59) 

(60) 

(61) 

(62) 

(6  V 

(64) 


Sib 
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FIG.  6.  Constant  cross  section  ion  distribution 
with  cold  neutrals. 


o-yro*)  «/0(CP)  +  AU*/,«r ), 

(65) 

which  yields 

A  U*/kTt=a/PkT.. 

Equation  (62)  is  integrated  to 

(66) 

"i=f  Mc)dv=  fl  /l+  P 

J-«  a  \  A/  an,  , 

) 

(67) 

w 

and  applied  to  the  Poisson  equation  (8)  to  produce 

,,+£H 

(68) 

Under  the  assumption  that  the  Debye  length  is  short  compared  to  the  ion  mean  free  path, 
(4i Ttfno/kT,  )2(  1/on.  )2»  1, 


Eq.  (68)  results  in 

P/<mn  » a/n.akT, (2  +  a/nmakT, )  (69) 

and 

W/kT,  -  1/(2  +  a/nmakT.).  (70) 

The  Bohm  criterion  is  satisfied  at  AU*  to  the  first  order  by  virtue  of  (68 ).  And  interestingly,  the  presheath  potential  rise  for 
a  *  0  is  exactly  that  required  by  the  cold  ion  Bohm  criterion.  Figure  6  presents  the  results  for  cold  neutrals  with  a/n„akT, 
»  I.  From  the  ion  distribution  at  AU*,  the  mean  ion  velocity  into  the  sheath  can  be  determined  to  be  v  *  1.0 6yjkT,/m, , 
while  the  e*act  solution  of  Rietnann  gives  C  «  1.27 yjkT,/mt-,  thus  the  first-order  asymptotic  result  appears  close. 

V.  CONCLUSIONS 

It  has  been  shown  that  approximate  collisions!  presheath  solutions  can  be  obtained  for  a  variety  of  collision  terms.  In 
particular  the  constant  collision  frequency  case  has  been  solved  approximately,  whereas  previous  attempts  at  exact  solutions 
have  found  this  case  intractable.  In  addition,  it  has  been  shown  that  higher-order  corrections  can  be  made  a  regular  and 
tractable  fashion.  Also  the  return  of  ions  from  the  collisionless  sheath  can  be  treated. 
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APPENDIX:  MULTIEXPONENTIAL  FORMULATION 

In  the  previous  sections  we  have  calculated  only  the  first-order  terms  in  the  ion  distribution  and  presheath  potential  rise. 
Also,  we  have  implicitly  made  the  same  first-order  approximation  for  electrons: 
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*#  *"o(!  —  AU/kT,).  (Al) 

A  complete  multiexponential  expansion  can  also  be  constructed  that  correctly  calculates  the  second-  and  higher-order  terms. 
Potential  in  the  presheath  is 

V  =  U0  +  AU  +  V 2  +  tij  A  U 3  +  •  •  • , 

where  U0=*ax and  AU »  txp(0x). Thus 


~^a+0AU  -I-  Ifia^AU*  +  30aAUi  + 
dx 


and 


d(AU) 


0AU 


dU  a  +  0AU  +  20a7AU2  -f  'i0a3AUi  -f 

which  transforms  the  Boltzmann  equation 

«L(vJL(v)!*£±±V(v))  =  m 

dxl  dAU  dU  X  m  dv  )  V*  )t 


into 


or 


±—  (a  +  0AU  +  WoAV7  + 
a  all  m 


l:  ± 


£'■>-(£)• 


iMf).],’ 

At/:  v0Uo)±L*k{v)±2.§L{v)  =  \(y.\  |  , 

m  m  mm;  lAdf  /day 

At/J:  2t;iS/J(l;)±^%l;)±^^-(l;)±i^(t,)«=[(f)  ]  , 

m  dv  m  dv  m  dv  Wdt  Je\Av' 


(A2) 

(A3) 

(A4) 

(A5) 

(A6) 

(A7a) 

(A7b) 

(A7c) 


At/-:  . 

m  dv  m  dv  m  dv  m  dv  l\«?r  JcUw 


(A7d) 


The  Poisson  equation  (8)  becomes 
0 zAt/+  (20)3«,At/2  +  (30)2a,At/3  +  ••• 

»4m^At/^J  /,  (v)dv  -  J  /„  (v)dvj 
+«c  fa(v)dv  —  J  +  "-j. 

(A8) 
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THE  TEC  PROGRAM  RESULTS 

The  TEC  program  results  shown  here  incorporate  the  asymptotic  presheath  work  and 
give  good  agreement  except  at  low  current  density.  This  disagreement  is  still  not  under¬ 
stood. 


-39- 


Comparison  of  TEC  Program  Results  to  Experimental  Data 

Reference:  Advanced  Thermionic  Technology  Program,  Vol.  4,  Oct.  1984.  Thermo 
Electron  Corporation,  DOE/ET/1 1292-2 


Output  Voltage  (volts) 
6^ 


Voltage 


TEC  INITIAL  DATA  SUMMARY 


PHYSICAL  OPERATING  CONDITIONS- 


1.50 


EMITTER  TEMPERATURE 

(TE)  = 

1700.0 

KELVIN 

COLLECTOR  TEMPERATURE 

(TC)  * 

773.0 

KELVIN 

EMITTER  WORK  FUNCTION 

(EWF)  = 

2.642 

EV 

COLLECTOR  WORK  FUNCTION 

(CWF)  * 

1.630 

EV 

CONVERTOR  PRESSURE 

(PN)« 

1.541 

TORR 

GAP  THICKNESS 

(D)« 

0.254 

MM 

OPERATING  CURRENT 

(J)» 

2.000 

AMPS/CMA2 

TEC  FUNCTION  SETTINGS 

DIAGNOSTIC  LEVEL 

(CHKDOT) 

«  1 

RESTART  SEQUENCE 

(OFILE) 

-  0 

POINT  DENSITY 

(N) 

*  11 

PHYSICAL  PARAMETERS  EVALUATED- 


RI CHARD SON  CURRENT  (JRIC) 

REFERENCE  DENSITY  (NR) 

CHARACTERISTIC  TIME  (TCHAR) 
NONDIM  CURRENT  (I) 

NONDIM  EMISSION  (ENR) 
KNUDSEN  NUMBER  (KN) 

SORT (MASS  RATIO)  (SMR) 
MEAN  FREE  PATH  RATIO  (LAMDAR) 


0.51E+01  AMPS/CMA2 
0.10E+15  1/CMA3 
0.0208  SECS*E-06 
0.0196 

0.008 (NRIC/NR) 
0.0791 
0.0020 
0.3344 


1.00 


0.50 


0.00 


0.00 


100.00  200.00  300.00  400.00 

T  in  Charcteristic  Time 

3/ 


500.00 


|0;u$  tc 


TEC  INITIAL  DATA  SUMMARY 


PHYSICAL  OPERATING  CONDITIONS 


EMITTER  TEMPERATURE 

(TE)  = 

1700.0 

KELVIN 

COLLECTOR  TEMPERATURE 

(TC)  = 

773.0 

KELVIN 

EMITTER  WORK  FUNCTION 

(EWF) = 

2.642 

EV 

COLLECTOR  WORK  FUNCTION 

(CWF) = 

1.630 

EV 

CONVERTOR  PRESSURE 

(PN)  = 

1.541 

TORR 

GAP  THICKNESS 

(D)  = 

0.254 

MM 

OPERATING  CURRENT 

(J)  = 

2.000 

AMPS /CM* 2 

TEC  FUNCTION  SETTINGS - 

DIAGNOSTIC  LEVEL 

(CHKDOT) 

=  1 

RESTART  SEQUENCE 

(OFILE) 

=  0 

POINT  DENSITY 

(N) 

*  11 

PHYSICAL  PARAMETERS  EVALUATED - 

■- 

RICHARDSON  CURRENT 

( JRIC) =  0. 

51E+01 

AMPS/CM*2 

REFERENCE  DENSITY 

(NR)*  0. 

10E+15 

1/CMA3 

CHARACTERISTIC  TIME  (TCHAR) = 

0.0208 

SECS*E-06 

NONDIM  CURRENT 

(I)* 

0.0196 

NONDIM  EMISSION 

(ENR) = 

0.008 (NRIC/NR) 

KNUDSEN  NUMBER 

(KN)* 

0.0791 

SQRT (MASS  RATIO) 

(SMR)  = 

0.0020 

MEAN  FREE  PATH  RATIO  (LAMDAR)  =  0.3344 

TIME  SETTINGS - 

NSTEPS*  1 

T2*  500.0 
DELTAT*  1.000 
DTP*  1.000 
LSF*  10 


RESULTS  AT  TIME  *  0.00 


ENE  - 

0.027 

ECHI* 

CNE  * 

0.000 

CCHI* 

PHIB* 

0.000 

VD  * 

# 

NDOT (#) 

NEB  (♦) 

0 

0.0015 

-0.154 

1 

-0.0151 

0.298 

2 

-0.0158 

0.661 

3 

-0.0141 

0.943 

4 

-0.0122 

1.145 

5 

-0.0104 

1.267 

6 

-0.0088 

1.309 

7 

-0.0072 

1.271 

8 

-0.0058 

1.153 

9 

-0.0050 

0.953 

10 

-0.0083 

0.670 

11 

-0.0340 

0.270 

OPERATING  VOLTAGE*  1.074 

.241  EALPHA*  0.460 
.734  C ALPHA*  0.482 
.422  EMISS  =.268E-0l 


TDOT (#) 

TAU ( # ) 

-0.0526 

1.52 

-0.0537 

1.52 

-0.0555 

1.51 

-0.0569 

1.51 

-0.0582 

1.50 

-0.0594 

1.50 

-0.0608 

1.50 

-0.0623 

1.49 

-0.0641 

1.49 

-0.0665 

1.48 

-0.0698 

1.47 

-0.0768 

1.45 

32 


12 

-0.0704 

-0.551 

-0.0807 

1.43 

RESULTS  AT  TIME  - 

100.00 

OPERATING  VOLTAGE-  0.982 

ENE  - 

0.098 

ECHI- 

3.183 

EALPHA- 

0.617 

CNE  « 

0.000 

CCHI- 

3.069 

CALPHA- 

0.719 

PHIB* 

0.000 

VD  =  • 

•0.203 

EMISS  = 

. 981E-01 

t 

NDOT(f) 

NEB  (#) 

TDOT(i) 

TAU  ( # ) 

0 

0.0003 

-0.035 

0.0001 

1.14 

1 

-0.0007 

0.081 

0.0001 

1.13 

2 

-0.0018 

0.200 

0.0000 

1.11 

3 

-0.0028 

0.316 

-0.0001 

1.10 

4 

-0.0037 

0.419 

-0.0001 

1.09 

5 

-0.0043 

0.501 

-0.0002 

1.08 

€ 

-0.0047 

0.552 

-0.0003 

1.07 

7 

-0.0046 

0.562 

-0.0003 

1.06 

8 

-0.0043 

0.524 

-0.0004 

1.05 

9 

-0.0035 

0.431 

-0.0005 

1.03 

10 

-0.0022 

0.281 

-0.0006 

1.01 

11 

-0.0006 

0.074 

-0.0008 

0.97 

12 

0.0013 

-0.170 

-0.0009 

0.95 

RESULTS  AT  TIME  « 

200.00 

OPERATING  VOLTAGE-  0.946 

ENE  - 

0.232 

ECHI- 

2.324 

EALPHA- 

0.587 

CNE  - 

0.000 

CCHI- 

2.167 

CALPHA- 

0.750 

PHIB- 

0.000 

VD  «  - 

•0.448 

EMISS  - 

.232E+00 

« 

NDOT ( # ) 

NEB (f ) 

TDOT ( # ) 

TAU  ( # ) 

0 

0.0001 

-0.014 

0.0006 

1.19 

1 

-0.0003 

0.034 

0.0008 

1.19 

2 

-0.0007 

0.084 

0.0008 

1.17 

3 

-0.0011 

0.134 

0.0008 

1.15 

4 

-0.0015 

0.179 

0.0007 

1.13 

5 

-0.0018 

0.216 

0.0006 

1.11 

6 

-0.0020 

0.240 

0.0005 

1.09 

7 

-0.0020 

0.247 

0.0005 

1.08 

8 

-0.0019 

0.232 

0.0004 

1.06 

9 

-0.0016 

0.193 

0.0002 

1.03 

10 

-0.0010 

0.127 

0.0001 

1.00 

11 

-0.0003 

0.033 

-0.0002 

0.93 

12 

0.0006 

-0.081 

-0.0003 

0.90 

RESULTS  AT  TIME  « 

300.00 

OPERATING  VOLTAGE-  0.910 

ENE  « 

0.526 

ECHI- 

1.440 

EALPHA- 

0.551 

CNE  - 

0.000 

CCHI- 

1.374 

CALPHA* 

0.754 

PHIB- 

0.000 

VD  -  • 

-0.698 

EMISS  - 

. 526E+00 

# 

NDOT (#) 

NEB  (#) 

TDOT  (#) 

TAU  ( # ) 

0 
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APPENDIX  A  -  Fokker  Planck  Collisions  Presheath 


This  theory  has  been  developed  under  this  grant  and  is  found  to  be  applicable  to  fully 
ionized  plasmas  but  was  not  incorporated  into  the  Thermionic  Convertor  work  due  to  its 
computational  complexity. 
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SUMMARY 

The  Maxwellian  sources  and  charge  exchange  terms  used  to  model  particle  inter¬ 
actions  in  current  presheath  models  do  not  represent  the  Coulomb  collisions  taking 
place  in  fully  ionized  plasmas.  These  models  approximate  the  collisional  effects  in 
the  presheaths  of  partially  ionized  plasmas  but  are  used  to  implicitly  extrapolate 
the  interesting  case  of  fully  ionized  plasmas.  The  present  study  uses  a  Fokker  - 
Planck  collision  term  which  models  the  limit  of  the  small  angle  Coulomb  collisions 
that  occur  in  fully  ionized  plasmas.  Normally  these  small  angle  collisions  dominate 
the  particle  interactions  of  fully  ionized  plasmas.  The  Boltzmann  equation  coupled 
with  the  Fokker  -  Planck  term,  and  the  Poisson  equation  have  been  expanded  using 
an  exponential  asymptotic  technique.  These  equations  have  been  solved  numerically 
to  determine  the  time  dependent  evolution  of  the  presheath.  The  results  presented 
show  the  presheath  potential  structure  and  particle  distribution  in  velocity  space. 
The  model  produces  a  self-consistent  and  accurate  potential  structure.  The  particle 
velocity  distribution  in  the  presheath  has  the  correct  acceleration  of  ions  toward  the 
wall  but  because  the  Fokker  -  Planck  collision  term  only  models  the  limit  of  small 
angle  collisions  it  is  unable  to  clear  the  particle  distribution  of  returning  ions.  The 
collisional  processes  become  dominated  by  the  effects  of  the  large  angle  collisions  as 
the  Debye  sheath  edge  is  approached.  This  study  has  found  that  a  presheath  model 
which  describes  the  Coulomb  collisions  occurring  in  a  fully  ionized  plasma  must 
account  for  both  the  small  angle  and  the  large  angle  particle  collisions  to  explain 
the  clearing  out  of  returning  ions  that  must  exist  for  the  transition  to  an  absorbing 
wall. 
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CHAPTER  I 
INTRODUCTION 

The  interaction  of  man  and  plasma,  in  some  form,  exists  at  almost  all  levels 
of  society.  A  plasma  is  an  ionized  gas  that  has  a  collective  behavior  in  an  electro¬ 
magnetic  field.  Plasmas  exist  in  everyday  devices  like  flourescent  lights,  neon  signs, 
and  electric  arc  welders.  An  understanding  of  the  basic  behavior  and  interaction  of 
plasmas  is  essential  to  the  advancement  of  all  current  plasma  applications  and  to 
the  discovery  of  new  applications.  This  thesis  involves  the  study  of  how  a  plasma 
interacts  with  the  walls  and  surfaces  with  which  it  comes  in  contact. 

Why  is  it  important  to  understand  plasma  -  wall  interactions?  Two  basic  reasons 
answer  this  question.  First,  a  plasma  has  a  strong  effect  on  any  surfaces  it  comes 
in  contact  with.  The  high  temperature  plasma  can  erode  or  destroy  any  surface 
quickly  pitting  and  changing  a  wall  which  may  need  to  maintain  a  particular  profile 
or  surface  condition.  Secondly,  the  wall  affects  the  characteristics  of  the  plasma. 
A  surface  r  an  have  a  profound  effect  on  the  plasma  depending  on  the  amount  and 
rate  at  which  it  can  absorb  energy.  Examples  of  situations  in  which  plasma  -  wall 
interactions  are  of  importance  include: 

•  Diverter  plates  in  magnetic  confinement  fusion  reactors. 

•  The  rails  in  a  plasma  rail  gun. 

•  Any  body  ( like  the  space  shuttle  )  upon  reentry  to  the  atmosphere. 

•  Plasma  switches. 


Plasma  etching. 
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•  Almost  every  other  use  or  occurrence  of  plasmas. 

This  study  Is  primarily  applicable  to  fully  ionized  plasmas.  The  hot  tempera¬ 
tures  necessary  to  produce  fully  ionized  plasmas  occur  only  in  situations  like  on  the 
surfaces  of  diverter  plates  in  Tokamak  fusion  reactors. 

The  development  of  a  mathematical  model  to  represent  the  plasma  -  wall  in¬ 
teraction  region,  or  sheath,  and  a  numerical  solution  to  this  model  is  the  focus  of 
this  thesis.  An  understanding  of  the  interaction  between  the  plasma  and  the  wall 
is  achieved  with  a  time  dependent  solution  to  the  sheath  region.  If  the  potential 
structure  and  the  particle  velocity  distributions  are  known  for  every  location  in  the 
sheath  then  the  energy  going  into  the  surface  can  be  determined.  In  this  way  the 
results  of  this  study  can  be  used  as  a  boundary  condition  for  problems  involving 
plasma  characteristics  and  for  problems  invloving  the  surface  physics  of  plasma 
devices. 
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CHAPTER  H 
BACKGROUND 

A  plasma  will  naturally  maintain  itself  in  a  neutral  and  field  free  state.  Ap¬ 
plication  of  forces  and  processes  that  try  to  alter  the  equilibrium  are  resisted  by 
the  plasma.  A  surface  within  a  plasma  that  is  not  at  the  same  potential  as  the 
plasma  will  be  shielded  from  the  remainder  of  the  plasma  by  a  sheath.  The  outer 
edge  of  this  sheath  is  nearly  at  the  plasma  potential.  Bohm11*  first  came  up  with  a 
criterion  to  determine  the  extent  of  the  sheath.  Bohm  modeled  the  sheath  region 
as  completely  collisionless.  He  also  considered  that  the  transition  region  from  this 
collisionless  sheath  to  the  plasma  was  too  small  to  be  important. 

More  recent  work  has  been  done  to  describe  this  transition,  or  presheath,  region. 
Self!3!  has  an  exact  solution  to  the  sheath  equation  and  has  shown  that  the  collision¬ 
less  sheath  makes  a  transition  directly  to  the  neutral  plasma  in  the  limit  as  ^  —  0, 
where  A o  is  the  Debye  length  and  L  is  the  plasma  dimension.  Emmert  et  aU3!  has 
determined  a  presheath  structure  based  on  the  assumption  of  a  Maxwellian  source 
of  ions  to  model  the  particle  collisions.  The  solution  to  this  model  shows  that  the 
transition  point  from  the  sheath  to  the  presheath  has  a  finite  electric  field  strength. 
Bissell  and  Johnson!*!  have  perfomed  a  similar  solution  using  a  Maxwellian  source  of 
ions.  In  contrast  to  Emmert  et  al.,  Bissell  and  Johnson  have  found  that  the  electric 
field  strength  becomes  infinite  at  the  sheath  edge.  This  solution  agrees  with  the 
fluid  and  cold  ion  models.  In  a  recent  paper  Bissell!8!  shows  that  Emmert  obtained 
a  finite  electric  field  strength  because  the  Maxwellian  source  term  used  produced 
no  ions  at  the  point  of  zero  velocity.  Bissell  and  Johnson  used  a  more  realistic 
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Maxwellian  source  that  produced  ions  at  the  zero  point  in  velocity  space  for  their 
solution. 

Another  approach  to  the  problem  involves  the  use  of  a  charge  exchange  term  to 
model  the  particle  collisions.  Riemannl6!  has  produced  results  using  this  technique. 
In  a  recent  paper  by  Main!7!  a  charge  exchange  model  is  used  to  obtain  a  solution 
to  the  presheath  potential  structure  and  particle  distribution.  This  model  involves 
an  asymptotic  approximation  of  the  plasma  equations.  The  Boltzmann  and  Poisson 
equations  are  asymptotically  expanded  and  then  solved  analytically  when  combined 
with  a  charge  exchange  model  of  the  particle  collisions. 

All  of  these  sheath  and  presheath  solutions  have  modeled  particle  collisions  by 
large  instantaneous  changes  in  particle  velocity.  These  models  do  not  represent  the 
Coulomb  collisions  occurring  in  the  presheath  of  a  fully  ionized  plasma. 

The  current  study  extends  the  asymptotic  solution  presented  by  Mainl7!  to  in¬ 
clude  a  Fokker  -  Planck  collision  term  instead  of  the  charge  exchange  term.  Unlike 
the  previous  collision  terms  used,  the  Fokker  -  Planck  term  describes  the  Coulomb 
collisions  that  exist  within  a  fully  ionized  plasma.  The  addition  of  the  Fokker  - 
Planck  term  necessitates  the  use  of  numerical  techniques,  rather  than  analytical 
techniques,  to  obtain  a  solution.  In  using  the  Fokker  -  Planck  term  the  collision  pro¬ 
cesses  are  being  modeled  directly.  The  model  developed  obtains  the  time  dependent 
evolution  of  the  presheath  for  a  fully  ionized  plasma. 
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CHAPTER  HI 

MODEL  FORMULATION 

3.1  Concepts 

In  order  to  have  a  complete  understanding  of  the  problem  at  hand  certain 
concepts  need  to  be  presented  which  will  help  in  understanding  the  overall  structure 
of  the  model. 

1)  Debye  Length  (A^)  -  The  shielding  distance  beyond  which  the  particle  charge 
effect  is  weak.  This  is  the  natural  charge  separation  distance.  Negatively  charged 
particles  become  surrounded  by  positively  charged  particles  and  vice  versa,  thus, 
balancing  the  overall  charge  at  any  point  (  see  figure  3.1 ).  There  is  a  point  beyond 
which  a  particle  is  not  effected  by  the  specific  charge  but  responds  to  the  influence 
of  the  entire  plasma.  The  thermal  effects  in  the  plasma  become  dominant  over  the 
electric  field  strength. 

2)  Mean  Free  Path  (A<)  -  The  average  distance  a  particle  travels  before  its 
trajectory  has  been  altered  by  ninety  degrees.  The  mean  free  path  is  a  function  of 
the  density  of  the  plasma.  The  denser  the  plasma  the  shorter  the  mean  free  path. 
For  the  plasma  under  consideration  in  this  study  A <  »  XD. 

3)  The  coordinate  system  used  to  describe  the  plasma  -  The  coordinate  system 
used  in  the  model  is  known  phase  space.  In  this  system  any  point  is  described 
using  three  position  coordinates  and  three  velocity  coordinates.  Any  orthogonal 
coordinate  system,  cartesian,  cylindrical,  spherical,  can  be  used  to  describe  both 
the  position  and  the  velocity  components. 
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4)  Distributions  and  Distribution  Functions  -  Plasmas  are  studied  in  a  collective 
sense.  The  motion  of  the  entire  plasma  and  not  individual  particles  is  described 
by  the  model.  Therefore,  the  velocity  of  the  plasma  at  any  given  location  must  be 
described  by  a  distribution.  The  distibution  function  describes  the  overall  particle 
velocity  distibution. 

5)  Potential  -  In  a  plasma  the  wall  potential  is  greater  than  the  neutral  plasma 
potential.  The  lighter,  thus,  faster  electrons  are  absorbed  by  the  wall  faster  than 
the  heavier  and  slower  ions.  A  net  positive  charge  exists  near  the  wall,  increasing 
the  potential  (  see  figure  3.2  ).  The  potential  at  the  physical  interface  between 
the  wall  and  the  plasma  is  dependent  on  the  rate  at  which  ions  are  absorbed  by 
the  surface.  In  this  study  u  =  -e<&  where  e  is  the  electron  charge  and  4  is  electric 
potential  in  electron  volts  so  that  U  has  units  of  energy.  The  addition  of  the  negative 
sign  defines  potential  in  the  reverse  of  the  usual  sign  convention  so  that  increasing 
potential  repels  electrons. 

6)  Collision  Possibilities  using  the  Fokker  -  Planck  Collision  term  -  To  describe 
the  overall  sturcture  of  the  sheath  the  various  collision  possibilities  must  be  included 
in  a  comprehensive  model.  The  Fokker  -  Planck  term  describes  the  four  major 
collision  possibilities. 

1)  Ion  -  Ion 

2)  Ion  -  Electron 

3)  Electron  -  Ion 

4)  Electron  -  Electron 

The  collision  model  does  not  take  into  account  thr  e  body  collisions.  Three 
body  collisions  are  very  rare,  as  such,  the  model  is  not  hampered  by  the  lack  of 
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terms  to  describe  these  collisions. 

3.2  Wall  Region  Model 

The  model  of  the  Plasma  -  Wall  region  can  be  broken  into  three  areas. 

1)  Neutral  Plasma  Region  ( 0{L ))  -  The  neutral  plasma  region  represents  the 
majority  of  the  system  and  can  be  considered  to  have  a  physical  width  that  is  on 
the  order  of  the  overall  dimension  of  the  system,  L.  This  region  is  considered  to  be 
fully  collisional.  The  velocity  distribution  is  near  Maxwellian  and  as  such  can  be 
modeled  by  fluid  type  equations  (  see  figure  3.3  ). 

2)  Debye  Sheath  Region  (O(Ap))  -  This  region  is  a  very  thin  area  directly  adjacent 
to  the  wall.  Its  width  is  considered  to  be  on  the  order  of  a  Debye  length  and  since 
A »  »  \D  no  collisions  are  expected  in  this  region.  This  collisionless  sheath  was 
first  modeled  by  Langmuir  1*1  and  Bohml1!  and  is  considered  very  well  known  and 
understood. 

3)  Collisional  Presheath  Region  (0(A,-))  -  This  is  a  transition  region  between  the 
collisional  neutral  plasma  and  the  collisionless  Debye  sheath  region.  It  is  considered 
to  have  a  physical  width  on  the  order  of  a  mean  free  path.  Therefore,  collisions  are 
expected  but  at  the  same  time  the  region  cannot  be  considered  fully  collisional. 

The  potential  must  transition  from  a  lower  level  in  the  neutral  plasma  to  a 
higher  level  at  the  wall.  The  goal  of  this  study  has  been  to  ohtain  a  time  dependent 
model  of  the  evolution  of  the  presheath  region  which  asymptotically  approaches  the 
known  potential  in  both  the  neutral  plasma  and  in  the  Debye  sheath  region. 

In  order  to  show  the  validity  of  the  three  region  model  an  example  of  Debye 
sheath  width  in  relation  to  the  overall  wall  region  is  appropriate.  For  this  example 
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average  hydrogen  fusion  plasma  characteristics  have  been  assumed: 


Ti  >  105K 


no  =  102Om~s 


If  the  Debye  length  within  this  plasma  is  calculated  an  order  of  magnitude  value 
for  the  Debye  sheath  width  is  determined.  An  appropriate  equation  for  the  Debye 


length  in  meters  is:M 


From  this  equation: 


XD  =  2.18  X  10-8m 


The  overall  sheath  width  is  on  the  order  of  a  mean  free  path.  An  appropriate 
equation  for  the  mean  free  path  in  meters  is:l10l 

For  a  singly  ionized  plasma  Z  »  1.  Using  this  equation  and  the  above  example 
plasma  characteristics  the  mean  free  path  can  be  calculated. 


A,-  -  1.14  x  10~*m 

This  is  an  order  of  magnitude  estimate  value  for  the  width  of  the  entire  sheath 
region.  Since  the  Debye  sheath  width  is  on  the  order  of  a  Debye  length  it  can  be 
seen  that  the  collisionless  sheath  is  very  thin  in  comparison  with  the  entire  wall 

region. 

In  order  to  obtain  an  idea  of  the  importance  of  the  electric  field  in  the  wall 
region  an  order  of  magnitude  analysis  is  useful.  The  magnitude  of  the  electric  field 
is  proportional  to  the  thermal  energy  per  length  scale. 

E  ~ 

z 


(3.3) 


In  the  sheath  region  the  length  scale  is  the  Debye  length. 


E  ~ 


kT\ 

Xd 


(3.4) 


Therefore,  in  this  region  the  electric  field  is  very  significant  since  XD  is  very  small. 
The  collisional  effects  are  small  in  comparison,  and  can  be  neglected. 

In  the  presheath  region  the  length  scale  is  the  mean  free  path. 


E  ~ 


Xi 


(3.5) 


Therefore,  the  electric  field  strength  is  on  the  order  of  the  collisional  effects  making 
both  important  factors  within  this  region. 

In  the  neutral  plasma  region  the  length  scale  is  the  overall  system  dimension. 

(3.6) 


Therefore,  the  electric  field  is  very  weak  and  can  be  neglected  in  comparison  with 
the  collisional  effects. 


&3.-Pr£3heath  Modal 

3.3.1  Equations  Describing  the  Collective  Behavior  of  a  Plasma 

The  primary  equation  used  to  describe  the  behavior  of  a  plasma  is  the  Boltz¬ 
mann  equation.  The  Boltzmann  equation  represents  the  collective  motion  of  many 
charged  particles  moving  in  an  electromagnetic  fieldln>. 


dt  dx  m  dx  9v  '  dt '  • 


(3.7) 


Where  the  +  sign  is  for  ion  particles  and  the  -  sign  is  for  electrons.  In  the  Boltzmann 
equation  /  is  the  particle  distribution  function,  and  is  defined  such  that  n  -  fdv. 


Tire  quantity  t  is  time,  v  is  a  velocity  vector,  x  is  a  position  vector,  m  is  the  particle 
mass,  and  V  is  the  potential.  The  term  on  the  right  hand  side  represents  particle 
collisions  and  can,  take  on  various  forms  depending  on  the  model  being  used. 

Another  important  equation  in  describing  the  bahavior  of  a  plasma  is  the  Pois¬ 
son  equation.  The  Poisson  equation  is  the  elementary  definition  of  potential  as  the 
collective  effect  of  charged  particles  on  a  point!7!. 

q3^J  fi{v,AU)dv- J  /,(v,AV)<ivj  (3.8) 

Where  q  is  the  elementary  charge,  the  subscript  i  refers  to  ions  and  the  subscript 
*  refers  to  electrons.  The  term  in  brackets  is  the  ion  -  electron  density  difference. 
The  potential  is  the  driving  force  in  the  Boltzmann  equation.  The  Poisson  equation- 
relates  the  potential  to  the  particle  distribution. 

To  complete  the  set  of  equations  necessary  for  a  full  description  of  the  presheath 
a  collision  term  must  be  chosen  to  model  the  particle  interactions.  This  study  uses 
the  Fokker  -  Planck  collision  term  to  model  the  particle  collisions.  The  Fokker  - 
Planck  term  represents  the  right  hand  side  of  the  Boltzmann  equation!11!. 


-  rf  d  (f  m  3  ^j)  i  1  32  (f  d*g  >1 
3t  Je  dvi^  2M  dvi  '  2dvidvk''/  dvjdv),* } 


Where, 


r  =  ?2|?rtInA 
4xt%m3 

9(v)  *  j  fW)  (  v  -  v'  |  dvf 


m  +  mf 


(3.9) 


(3.10) 

(3.11) 

(3.12) 


Where  M  is  called  the  reduced  mass.  No  superscript  refers  to  the  particle  species 
undergoing  the  collisions  and  the  superscript  /  refers  to  the  particle  species  with 
which  the  collision  occurs. 
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The  Fokker  -  Planck  equation  describes  the  Coulomb  collision  between  two 
charged  particles.  Certain  restrictions  and  assumptions  are  made  when  using  the 
Fokker  -  Planck  collision  term.  First,  it  best  describes  fully  ionized  plasmas.  The 
collision  term  models  charged  particle  interaction  and  is  most  accurate  for  plasmas 
with  few  neutrals.  This  situation  occurs  only  on  the  hottest  of  plasma  surfaces  like 
the  diverter  plates  in  Tokamak  fusion  reactors. 

The  second  restriction  involves  the  type  of  collisions  that  the  Fokker  -  Planck 
term  models.  The  overwhelming  majority  of  particle  collisions  lead  to  only  small 
deflections  in  the  particle  trajectories.  The  Fokker  -  Planck  term  describes  the  limit 
of  these  small  angle  deflections.  Finally,  the  model  does  not  take  into  account  three 
body,  and  higher  order,  collisions. 

3.3.2  Solution  Conditions 

The  three  equations  presented  in  the  previous  section  in  conjunction  with  the 
asymptotic  forms  of  potential  and  velocity  distribution  provide  the  necessary  infor¬ 
mation  to  determine  the  presheath  structure  if  two  additional  conditions  are  met. 

First,  if  the  equations  are  written  in  cylindrical  coordinates  the  particle  velocity 
distribution  is  axially  symmetric.  There  is  no  theta,  8 ,  dependence  of  the  velocity 
distribution.  Cylindrical  coordinates  are  used  for  both  the  velocity  and  the  position. 
The  V  direction  is  perpendicular  to  the  wall  (see  Figure  3.4)  with  the  positive 
direction  being  defined  into  the  wall.  The  coordinates  R,8 ,*  have  been  used  in 
velocity  space  for  convenience. 

The  second  condition  for  a  solution  to  these  equations  involves  an  assumption 
of  the  particle  velocity  distribution  parallel  to  the  surface.  For  this  model  the 
radial  velocity  distribution  has  been  assumed  to  take  the  form  of  a  Maxwellian 
distribution.  In  addition,  the  temperature  in  the  radial  direction  has  been  assumed 
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to  be  uniform  and  constant  at  all  V  locations.  Thus,  the  radial  velocity  distribution 
is  constant  for  any  position.  Figure  3.5  is  a  schematic  of  these  conditions.  Note 
that  at  any  V  location  and  rotation  angle,  $,  the  radial  velocity  distribution  is 
constant  and  follows  a  Maxwellian  distribution.  This  represents  the  conditions  of 
uniform  temperature  and  axial  symmetry,  throughout  the  wall  region.  Figure  3.5 
also  shows  a  representation  of  the  point  of  no  returning  ions.  This  is  the  point 
where  the  presheath  transitions  to  the  collisionless  sheath. 

The  conditions  of  uniform  temperature  and  radial  Maxwellian  distribution  al¬ 
though  good  approximations  are  not  exact  models  of  the  real  situation. 

The  overall  problem  reduces  to  one  dimension,  the  *  direction,  with  the  above 
conditions.  The  6  dependence  having  been  removed  by  the  axial  symmetry  and  the 
radial  dependence  having  been  removed  by  the  Maxwellian  assumption.  This  one 
dimensional  problem  can  be  solved  by  staight  forward  numerical  techniques. 

3.3.3  Expansion  of  the  Boltzmann  Equation 

The  presheath  model  involves  the  expansion  of  the  potential  and  velocity  into 
asymptotic  approximations.  The  potential  is  assumed  to  follow  an  exponential 
asymptotic  form. 

U  =  Uo  +  aiAU  +  ajAU7  +  •  •  •  (3.13) 

where 

U0  —  a*  and  A U » (3.14) 

«t»*s  are  parameters  which  describe  the  potential  structure.  Alpha,  a,  is  non-zero 
for  a  non-zero  potential  gradient  in  the  neutral  plasma.  A U  is  called  the  potential 
expansion  parameter. 

The  particle  distribution  in  velocity  space  is  a  function  of  potential  and  can  be 
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similarly  expanded. 


f(v,AU)  =*  /o(w)  +  AU  fx{v)  +  AU2  f7[v)  + 
The  Boltzmann  equation  can  be  written  as 


df  df  a[AU)  iaudf  .df* 

■  +  -  x  m  a,  a,,  t gt ** 


(3.15) 


(3.16) 


dt  '  d(AU )  dx  m  dx  dv 
where  ±  sign  is  respectively  positive  for  ion  particles  and  negative  for  electrons. 
The  potential  U  has  units  of  energy  and  is  defined  as  shown  in  figure  3.3  such  that 
U  »  — e*  where  4  is  electric  potential.  The  Boltzmann  equation  can  be  expanded 
using  equations  3.13  and  3.14.  In  addition  since  the  solution  is  one  dimensional  in 
velocity  space  the  velocity  derivatives  reflect  only  the  V  direction.  The  following 
expansions  are  used. 

<3I7) 

df  -  /i+2Atf/3+3Atfa/3  +  -.-  (3.18) 

PAU  (3.19) 


d(AU) 

d(AU) 


dx 


orr 

^  »  a  +  paxAU  +  20a,A U3  +  3pa3&U*  + 

^  +  AC r|£-  +  Atfa^  +  --- 
dx  ox  dx  dx 


(3.20) 

(3.21) 


Using  these  expansions  the  Boltzmann  equation  can  be  broken  down  by  order  in 
AU  assuming  the  collision  term  can  be  likewise  broken  down: 
dfo  a  dfo  f/£Z\  I 

dt  mdx  L1  at 


act : 

dt  P  11  m  dx  m  dx  [X  dt  '*J 

^  +  2 f.r,±z!h±hi!k.± 

dt  m  dx  m  dx 


(3.22a) 

[«'-L 

(3.22 4) 

2 Pa?  dfo 
m  dx 

(3.22c) 

AIT 


3f n  ,  _a_,  J.  a  dfn  ^  pax  dfn.x  ,  2pa7  Bfn-x 


+  nPxfn  ±  —  —  ± 
dt  max  m 


dx 


m  dx 


^  (n  -  l)gaw.t  dfx  nPdn  dfo 

m  dx  m  dx 


[<S>. 


(3.22d!) 


ay* 


sv 


14 


The  above  exponential  expansion  is  the  only  known  way  to  break  apart  the  Boltz¬ 
mann  equation.  To  complete  the  expansion  the  Poisson  and  Fokker  -  Planck  terms 
must  be  similarly  broken  down. 

At  this  point  it  is  appropriate  to  nondimensionalize  the  expansion  of  the  Boltz¬ 
mann  equation.  The  following  nondimensional  quantities  have  been  used. 


P.  a*  J 
'  TXR*3/* 

for  0  <;  t  <  n 

(3.23a) 

r 

ai/a E 

(3.236) 

z  =va7* 

(3.23c) 

aa 

(3.23d) 

mri^r 

**  nRr 

(3.23e) 

Pi  = 

mnK  V 

for  1  <  *  <  n 

(3.23/) 

Where  nA  is  a  reference  density,  a  and  a!  are  the  inverse  square  of  the  thermal 
velocities  of  the  colliding  particles,  a  *  and  o'  =  The  quantities  <p-\  through 
<pn  represent  the  potential  structure  of  the  presheath.  Fi  is  nondimens ionalized 
such  that  1  *  f?„FdZ.  Using  these  nondimensionalizations,  the  expansion  of  the 
Boltzmann  equations  becomes 


1: 
LU : 


At/ 


(3.24a) 

(3.246) 

(3.24c) 


dZ 


dZ 


,  ,  dFi  dFo  tBF\  1 


(3.24<f) 


S*Sf 
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3.3.4  Expansion  of  the  Poisson  Equation 


The  Poisson  equation, 

^-y=4?r q2^J  fi(v,AU)dv-  J  ft(v,AU)dv  (3.8) 

is  broken  down  using  the  same  technique  as  the  Boltzmann  equation.  Since: 

a3  tj 

^■4  -  02aiAU  +  402O2AU2 +902a$AU*  +  •••  (3.25) 

oxd 

[  f(v,AU)dv=*  f  fo(v)dv  +  AU  f  fi(v)dv  +  AU3  f  /a(«/)<f«H -  (3.26) 

J  —  00  7-00  /-oo  y-oo 

Using  these  expansions  the  Poisson  equation  can  be  broken  down  by  order  in  AU. 


1: 

fiO  dv-  J  /^(v)  <fuj 

(3.27a) 

AU: 

02ax 

f  fii[v)dv-  f  fti(v)dv 

-OP  J  —OP 

(3.276) 

AU1  : 

4/S2aa  =  4 xq2  ^ 

J  fit[v)dv-  j  /«a(w)rfvj 

(3.27c) 

AU"  :  n202an  =  4 xq2  ^  j  /,*(v)  dv-  J  fm(v)  <fwj  (3.27<f) 

The  assumption  of  Boltzmann  electrons  is  made  to  enable  the  numerical  calcula¬ 
tions  to  proceed  with  time  steps  on  the  order  of  an  ion  characteristic  time.  In  the 
asymptotic  presheath  the  Boltzmann  electron  assumption  becomes 

n.  =  noe-[««At/+**Atr*+-+-.AU‘)rf7  (3 .28) 


where  n,  is  the  electron  density,  T,  is  the  electron  temperature,  and  no  is  the  electron 
density  in  the  asymptotic  presheath  at  AU-O.  Expansion  of  (3.28)  in  terms  of  AU 
yields 
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With  the  assumption  that  \d  <  A*  the  Poisson  equation  reduces  to  equating  electron 
and  ion  densities  in  order  of  AU.  The  Poisson  and  Boltzmann  electron  equations  are 
nondimensionalized  in  the  same  manner  as  the  Boltzmann  equation.  The  nondi¬ 
mens  ionalized  Poisson  equation  (3.27a-d)  is  combined  with  the  Boltzmann  electron 
equation  (3.28)  to  become 


fm  F0dZ 


(3.30a) 

(3.306) 

(3.30c) 


where  »j  »  ^  and  may  be  specified  as  a  function  of  time.  This  equation  can  be  used 
to  solve  for  the  potential  structure  at  each  time  step.  2<  is  the  ion  temperature  and 
T,  is  the  electron  temperature. 

3.3.5  Expansion  of  the  Fokker  -  planch  Term 

The  Fokker  -  Planck  term  must  also  be  expanded  in  order  of  AU  but  first  must 
be  put  into  cylindrical  coordinates.  In  addition,  the  assumption  of  axial  symmetry 
must  be  accounted  for  ir.  the  term.  This  can  be  accomplished  by  expanding  each 
term  in  the  general  Fokker  -  Planck  term. 


«r 


(/— — v3*)  +  I_i!_ 

U  2 Mdvt  91  2  dvtdvu 


The  first  term  can  be  rewritten: 


(3.9) 


2  M 


wa  ^  9  (  dg\  1  d*g  d*g 

9  RdR\  dRJ  R*  d6*  dz* 


(3.31) 


t 

•it 


u 


(3.32) 
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(3.38) 
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With  this  definition  the  Fokker  -  Planck  term  becomes: 


Kir}.  -  i£r(/jr  &</<*+<)> ' « '  -O-ioa/l  IK'<‘+{»  i  ■ « i  ■ «)  <3m> 


For  brevity  let  /  denote  a  derivative  with  respect  to  t.  Using  this  notation  the  Fokker 
-  Planck  term  becomes: 


M30.  ■  'S('C  i « i  ■ «)  -  tA'&C  n-+() '  ■ 5 1 «)  <3-4°) 

The  Fokker  -  Planck  term  is  nondimensionalized  using  the  same  variables  as  the 
Boltzmann  equation. 

{If}. - ttMjLI nz+ou\ *) *Mr=S}£' r"iz+a in*)  o«) 

Let: 


+  (3  42) 

B(F)  /_”  r"(Z  +  f)  I  f  I *  (*•«) 

The  Fokker  -  Planck  term  can  be  written  in  a  compact  form  using  these  defined 
functions. 

{If  <s“> 

The  Fokker  -  Planck  term  can  be  expanded  by  order  in  LU  using  the  same  technique 
as  the  Boltzmann  equation. 

{If}.-  [w  (w.l)  +  £  (W.))] 

+  4I/[^i(F,X(Fo))  +  ^(f,B(fo))  +  Mr,))  +  ^(ftWl))] 

+  ••• 

+  ACT  [  £  +  £(*.-«*(*.)))] 


fa 


(3.45) 


19 


This  expansion  can  be  combined  directly  with  the  expansion  of  the  Boltzmann 
equation. 

3.3.6  Solution  Approach 

To  obtain  a  time  dependent  solution  to  the  Boltzmann  equation  with  the  Fokker 
-  Planck  collision  term  a  numerical  technique  is  necessary.  The  solution  presented 
has  been  limited  to  ion  -  ion  collisions  because  these  collisions  represent  the  majority 
of  the  coilisional  energy  and  momentum  transfer  within  the  presheath.  The  positive 
signs  in  the  Boltzmann  equation  must  be  applied  for  ion  -  ion  collisions.  The  ratio 
of  particle  mass  to  reduced  mass,  must  be  equal  to  two  for  like  particles.  The 
ratio  of  the  inverse  squares  of  the  particle  thermal  velocities,  must  be  one  for 
like  particle  collisions.  The  nondimensional  Boltzmann  equation  reduces  to  the 
following  form. 

(**<*>)  -  jz  (*B™) + *»•«•< +*>-fy 

*  ap  <  /  ,2  ,  .  a  ,  .  \  (3*^6) 

—  E  E 

imI  «n*l  '  ' 

where  the  summations  are  taken  to  be  zero  if  i  —  0.  The  functions  A(Fm)  and  B{Fm) 
can  be  written: 


1  -«• 


The  V  equations  in  the  expansion  are  solved  to  obtain  the  time  dependent 
particle  velocity  distribution.  The  Poisson  equation  is  employed  at  each  time  step 
to  obtain  the  potential  structure.  The  ratio  of  the  higher  order  equations  with 
respect  to  the  first  order  equation  eliminates  the  term  from  the 

Poisson  equation.  Using  this  technique  each  successive  component  of  the  potential 


20 


structure  can  be  determined  from  the  previous  components. 
XHafWg  ft  (ZZ\ 

‘/r»^o  dZ~  <Po\Tt) 

jr^F3dz 

ST^FodZ  <Po\T'J  2\<poJ  \TtJ 

I?,*  FjdZ  p*m\.  Vi£l(n±\a 

I-C  F0dZ  ~~  vo\T.  )  po<pq\T.J  G\<poJ 


(3.49a) 

(3.496) 

(3.49e) 


Using  these  equations  the  particle  velocity  distribution  and  the  potential  structure 
of  the  presheath  are  determined  as  a  function  of  time. 

Of  interest  in  this  study  is  the  point  at  which  there  are  no  returning  ions. 
This  is  the  presheath  -  sheath  interface.  This  occurs  when  the  net  ion  flux  away 
from  the  wall  is  zero.  To  calculate  this  point  in  the  presheath  it  is  necessary  to 
obtain  a  value  for  the  potential  expansion  parameter  L.U  such  that  when  the  overall 
particle  distribution  is  reconstructed  from  the  various  terms  in  the  expansion  no 
returning  ions  are  present.  Thus,  at  the  critical  point  of  no  returning  ions  the 
model  determines  the  total  particle  velocity  distribution,  the  potential  structure, 
and  a  value  for  the  potential  expansion  parameter.  From  the  potential  structure  and 
the  potential  expansion  parameter  the  presheath  height  at  the  point  of  no  returning 
ions  can  be  determined  (see  figure  3.5). 
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CHAPTER  IV 
NUMERICAL  1ECHNIQUE 
4.1  Problem  Approach 

The  solution  of  the  Boltzmann  equation,  as  written  in  equation  3.46,  coupled 
with  the  Poisson  equation  (3.49)  is  the  goal  of  the  numerical  procedure. 

The  general  approach  is  to  solve  the  Boltzmann  equation  for  the  particle  velocity 
distribution  using  a  partially  implicit,  partially  explicit  scheme.  Each  step  in  time 
the  Boltzmann  equation  is  solved  using  some  results  from  the  previous  time  step. 
In  equation  3.46  the  left  hand  side  is  solved  implicitly  while  the  right  hand  side  is 
solved  explicitly. 

77  -  sfi  (W.l)  -  £(w.>) +*>-1$ 
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The  left  hand  side  of  this  equation  can  be  put  in  a  matrix  form. 

(4.1) 

In  this  form  the  matrix  T(r)  is  an  m  x  m  matrix  created  from  the  left  hand  side  of  the 
Boltzmann  equation.  The  quantity  m  is  the  number  of  divisions  in  the  velocity  space 
2  chosen  for  the  numerical  scheme.  The  matrix  T{r)  is  computed  from  A(fo(r))  and 
B(Fo{r)).  The  values  of  these  derived  functions  are  taken  from  the  solution  to  the 
particle  velocity  distribution  at  the  previous  time  step,  r.  Numerical  derivatives  are 
used  to  represent  the  partial  derivatives  in  the  equation.  This  procedure  produces 
a  diagonal  matrix  where  all  elements  except  those  on  an  odd  number  of  centered 
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diagonals  are  zero.  The  number  of  diagonals  reflects  the  order  of  accuracy  in  the 
solution.  Diagonal  matrices  of  this  form  are  easily  and  quickly  inverted.  The 
Fi{r+ Ar)  matrix  is  an  mx  1  matrix  of  unknowns  that  represents  the  particle  velocity 
distribution  at  the  current  time  step.  V  equations  of  this  form  can  be  written 
corresponding  to  the  number  of  terms  in  the  expansion. 

The  right  hand  side  of  the  Boltzmann  equation  can  also  be  put  in  a  matrix 
form. 
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The  scalar  <pn  values  are  unknowns  and  represent  the  nondimensional  coefficients  in 
the  asymptotic  potential  structure  of  the  presheath. 

The  v?{t)  matrices  are  m  x  l  matrices  which  are  comprised  of  the  partial  deriva¬ 
tives  of  the  velocity  distribution  at  the  previous  time  step,  r.  They  represent  the 
first  summation  on  the  right  hand  side  of  the  Boltzmann  equation. 


m- 


az 


) 


The  dk  matrix  isanmxi  matrix  comprised  of  the  second  summation  on  the  right 
hand  side  of  the  Boltzmann  equation. 


ntal  '  ' 

All  values  of  the  distribution  and  the  functions  A(/m)  and  B[Fm)  are  taken  at  the 
previous  time  step,  r. 

Putting  together  equations  4.1  and  4.2  a  matrix  form  of  the  Boltzmann  equation 
is  created  that  can  be  solved  for  the  particle  velocity  distribution. 
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V  equations  of  this  form  can  be  written  corresponding  to  the  number  of  terms  in 
the  asymptoic  expansion  being  used. 

These  equations  are  quickly  inverted  to  obtain  the  particle  velocity  distribution 
at  the  current  time  step. 
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Where  the  Vn  and  D  matrices  represent  the  m  x  l  solution  matrix  to  the  inversion 
of  the  T(r)  matrix  with  the  corresponding  «„  or  d  matrix. 

The  particle  velocity  distribution  is  obtained  from  this  equation  using  the  <pn 
values  from  the  previous  time  step. 

Equation  3.49  is  employed  to  obtain  the  values  at  the  current  time  step. 
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The  new  distributions  are  integrated  numerically  and  the  <px  through  <pn  scalars 
are  determined  consecutively.  The  value  of  <po  is  input  and  is  a  nondimensional 
representation  of  the  coefficient  0  in  the  potential  expansion  parameter,  AU . 

For  each  time  step,  the  overall  particle  velocity  distribution  can  be  determined 
at  any  location  from  the  original  expansion  once  it  has  been  nondimensionalized. 

F{Z,  AIT)  rn  Fo{Z)  +  AITFX{Z)  +  AlT3F7{Z)  +  •  •  •  (4.5) 
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Where  the  •  quantities  are  nondimensional.  x  has  been  nondimensionalized  with 
respect  to  an  ion  mean  free  path. 


n/*r 


The  potential  structure  of  the  presheath  is  detemined  from  the  nondimensional 
form  of  the  original  expansion  of  potential. 


ir  =  u;  +  px ait  +  ^  a  it*  +  ■■■ 


(4.7) 


Using  this  procedure  the  time  dependent  evolution  of  the  presheath  is  obtained. 

The  point  of  no  returning  ions  occurs  where  the  integral  of  the  left  half  plane 
of  the  total  particle  velocity  distribution  is  zero. 

0»  [  F(Z)dZ  (4.8) 
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This  equation  can  be  rewritten  using  the  expansion  of  the  particle  distribution. 

0-  f°  F0(Z)dZ  +  ACT  f°  Fi(Z)dZ  +  AIT7  [°  F2(Z)dZ  +  •••  (4.9) 

e-oo  j  —  OO  J  —  oo 

Equation  4.9  can  be  solved  for  A  IT  Since  the  particle  distributions  are  now  known 
as  a  function  of  time  and  velocity.  The  entire  solution  at  the  point  of  no  returning 
ions,  is  known  with  this  last  piece  of  information. 


4.2  Numerical  Integration,  Differentiation,  and  Matrix  Inversion 

In  order  to  obtain  a  solution  to  the  potential  structure  and  particle  distribution 
in  the  presheath  it  is  necessary  to  develope  the  applicable  mathematical  tools.  The 
primary  techniques  needed  are  integration,  differentiation,  and  matrix  inversion. 
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4.2.1  Numerical  Integration 

Throughout  the  solution  integration  is  computed  using  a  Simpson’s  £  rule 
technique!12!. 

J  /(*)dl  =  J  (/l  +  4/a  +  2/s  +  4/4  +  2/5  + - 1-  2fn-\  +  4/n  +  /„+ 1)  (4-10) 

Where  h  is  the  spacing  between  the  points  and  /1  through  fn+i  represent  the  function 
values  at  each  point.  This  procedure  has  a  global  error  of  0(h4).  If  the  step  size  is 
chosen  appropriately  this  procedure  is  very  accurate. 

4.2.2  Numerical  Differentiation 

The  technique  for  determining  numerical  differentiation  is  a  second  order  ac¬ 
curate  scheme.  This  reduces  the  number  of  computations  while  maintaining  high 
accuracy.  Second  order  accurate  numerical  differentiation  requires  that  only  three 
points  be  known.  Thus,  the  ’T(r)'  matrix  contains  only  three  diagonals.  If  third 
order  accuracy  was  used  the  ’T(r)’  matrix  would  require  five  diagonals  to  represent 
the  five  points  needed  for  the  differentiation.  In  addition,  to  maintain  uniformity  a 
central  difference  technique  is  desirable  on  as  many  points  as  possible.  The  greater 
the  number  of  points  needed  for  each  derivative  the  more  points  that  require  forward 
or  backward  difference  techniques  (  rather  than  the  central  difference  technique). 
Below  is  a  list  of  the  techniques  used  to  obtain  derivatives!12!. 

Central  Difference 


dF  =  f(«+l)-f(s-l) 
dx  2h 

d3F  =  F(z  -t- 1)  —  2F[x)  +  F{x  —  1) 
dx 3  “ 

a*F  _  F(x  +  2)  -  2 F(x  +  1)  +  2^(1  -  1)  -  F{x  -  2) 
dx3  ~  2 h3 


(4.11a) 

(4.116) 


(4.11c) 


Forward  Difference 


-F(x  +  2)  +  4 F[x  +  1)  -  3F(x) 

—  a 

dx 

2  h 

32F 

F [x  +  2)  -  2F{x  +  1)  +  F (x) 

dx2 

h2 

3*F 

F(x  +  3)  -  3 F(x  +  2)  +  3F(x  +  1)  -  F(x) 

3s3 

h3 

Backward  Difference 


dF  3 F(x)  -  4 F(x  -  1)  +  3 F[x  -  2) 
dx  ~  2h 

32F  _  F(x)  -  2F{x  -  1)  +  F[x  -  2) 
dx1  h 2 

F(x)  -  3 F[x  -  1)  +  3F(x  -  2)  -  F{x  -  3) 
dx3  k3 


(4.12a) 

(4.126) 

(4.12c) 


(4.13a) 

(4.136) 

(4.13c) 


Where  h  is  the  grid  spacing.  The  derivatives  are  being  taken  about  point  x. 

It  is  worth  noting  that  the  third  derivative  equations  require  up  to  five  points. 
There  is  no  second  order  accurate  numerical  third  derivative  representation.  These 
equations  are  third  order  accurate.  This  does  not  effect  the  ’r(r)’  matrix  in  that 
it  contains  no  third  derivatives.  The  solution  procedure  requires  third  derivatives 
only  in  the  determination  of  the  function  B(Fi{r)). 

These  equations  are  used  throughout  the  solution  for  derivatives  with  respect 
to  veJ.  /city,  Z,  and  time,  r. 


4.2,3  Matrix  Inversion 


In  order  to  obtain  a  solution  a  procedure  for  inverting  a  diagonal  matrix  is 
necessary.  The  procedure  used  will  invert  any  centered  diagonal  matrix.  For  the 
second  order  accurate  case  the  matrix  in  question  is  tridiagonal.  The  procedure 
uses  Guassian  elimination  on  all  terms  below  the  center  diagonal  and  then  through 
back  substitution  determines  the  solution  vector.  This  technique  can  quickly  invert 
a  200  x  200  tridiagonal  matrix. 


4.3  Obtaining  a  Solution 
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As  in  any  numerical  model  certain  restraints  and  conditions  must  be  met  to 
obtain  an  accurate  solution.  This  model  requires  some  form  of  input  distribution 
function  and  in  order  to  obtain  higher  order  terms  must  also  have  a  perturbation 
applied  to  the  potential  structure.  In  addition,  certain  numerical  techniques  have 
been  used  to  remove  instabilities  in  the  model. 

4.3.1  Initial  Distribution 

To  model  the  presheath  region  an  initial  particle  velocity  distribution  that  con¬ 
forms  to  a  Maxwellian  profile  has  been  used.  This  profile  represents  the  distribution 
that  naturally  occurs  in  the  neutral  plasma  region.  The  idea  is  that  the  time  de¬ 
pendent  evolution  of  the  distribution  will  change  from  a  Maxwellian  at  time  zero  to 
a  shifted  new  form  as  the  presheath  is  entered.  The  Maxwellian  profile  is  initially 
given  to  the  zero  order  term  having  set  the  initial  conditions  of  all  higher  order 
terms  to  zero. 

If  the  potential  structure  of  the  presheath  is  not  perturbed  in  some  manner  then 
the  model  represents  the  neutral  plasma  region  and  the  particle  velocity  distribution 
remains  Maxwellian  (  as  it  should  ).  If,  however,  a  small  perturbation  in  the 
potential  atucture  is  added  (ie.  a  nonzero  ai,a3  ••  •)  then  the  model  readjusts  to 
describe  the  presheath  region.  In  this  manner  the  model  is  used  to  give  the  time 
dependent  evolution  of  the  presheath. 

4.3.2  Instability  Damping 

By  the  nature  of  the  implicit  -  explicit  technique  being  employed  certain  numer¬ 
ical  problems  are  expected  to  appear.  This  model  is  no  exception.  Two  techniques 
have  been  used  to  remove  these  instabilities. 
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The  most  important  thing  to  do  to  avoid  numerical  problems  in  a  scheme  of  this 
nature  is  to  ensure  that  as  much  as  possible  of  the  solution  is  computed  implicitly. 
In  addition,  once  some  new  data  has  been  calculated  it  should  be  applied  to  any 
new  calculations  immediately. 

In  this  model  each  term  in  the  expansion  of  the  particle  distribution  function 
uses  the  new  data  already  determined  in  calculating  all  of  the  lower  order  terms. 
Once  F0  is  determined  that  information  is  used  in  calculating  Ft.  This  idea  is 
repeated  for  the  higher  order  terms. 

A  second  method  applied  to  the  model  to  eliminate  oscilliatory  instabilities  that 
start  on  a  very  small  scale  and  grow  is  the  application  of  a  very  weak  averaging 
scheme  to  the  particle  distribution  functions.  Each  point  in  the  distribution  is 
weakly  averaged  with  the  points  on  either  side. 

F(Z)  -  a025jrfr  +  *>  +  FW  +  0  0 2sFlz  -  *)  (4.14) 

This  technique,  although  necessary,  has  the  negative  effect  of  falsely  increasing  the 
energy  in  the  system  by  spreading  the  distribution  slightly  (see  figure  5.1).  The 
change  is  very  small  and  can  be  considered  insignificant  with  respect  to  the  overall 
solution. 


4.4  Program  Structure 

The  entire  program  has  been  written  in  FORTRAN  and  can  be  run  on  either 
an  IBM  PC  AT  or  on  the  CYBER  mainframe.  The  code  has  been  written  in  a  seg¬ 
mented  manner  that  easily  allows  one  section  to  be  altered  without  having  to  alter 
other  sections.  The  overall  structure  of  the  program  consists  of  three  initialization 
programs,  three  input  data  files,  the  main  program,  and  three  output  data  files. 
The  main  program  contains  a  driver  and  seventeen  subroutines.  Several  of  the  sub- 
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routines  perform  operations  that  are  used  throughout  the  main  code.  Figure  4.1  is  a 
diagram  of  the  structure  of  the  program.  The  flexibility  of  the  code  is  derived  from 
the  generalized  subroutine  structure  and  the  ability  to  enter  a  variety  of  input  vari¬ 
ables.  The  driver  keeps  track  of  time  and  maintains  the  overall  operating  structure 
of  the  solution  while  the  subroutines  perform  the  necessary  manipulations.  Below 
is  a  list  of  the  function  of  each  program,  data  file  and  subroutine. 

AVE  -  Subroutine  to  smooth  distributions  by  averaging. 

CONSERV  -  Subroutine  to  determine  conservation  of  energy,  momentum,  and 
particles. 

CONSOUT  -  Conservation  output  data  file. 

CRF  -  Particle  distribution  initialization  program. 

CRPHI  -  Potential  structure  initialization  program. 

DENSITY  -  Subroutine  to  solve  for  a  new  presheath  structure. 

FDATA  -  Initial  particle  distribution  data  file. 

FDl  -  Subroutine  to  find  first  derivatives. 

FD2  -  Subroutine  to  find  second  derivatives. 

FD3  -  Subroutine  to  find  third  derivatives. 

FINDA  -  Subroutine  to  determine  ’A*  function. 

FINDB  -  Subroutine  to  determine  ’5’  function. 

FPINIT  -  Primary  initialization  program. 

FPOUT  -  Output  particle  distribution  data  file. 

FPSHETH  -  Main  program  driver. 

GETAB  -  Subroutine  to  make  A  and  B  function  vectors. 

INITDAT  -  Initialization  data  file. 

MAKED  -  Subroutine  to  make  d  matrix. 
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MAKEF  -  Subroutine  to  read  initial  particle  distribution. 

MAKEPHI  -  Subroutine  to  read  initial  potential  structure. 

MAKET  -  Subroutine  to  make  T  matrix. 

MAKEV  -  Subroutine  to  make  v  matrix. 

MODIAG  -  Subroutine  to  invert  diagonal  matices. 

PHIDAT  -  Initial  potential  structure  data  file. 

PHIDOUT  -  Output  potential  structure  data  file. 

SIMPS  -  Subroutine  to  perform  Simpson’s  rule  integration. 

TOT  -  Subroutine  to  obtain  total  distribution  at  point  of  no  returning  ions. 
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CHAPTER  V 

RESULTS 

The  ratio  of  ion  temperature  to  electron  temperature,  £,  has  been  set  to  one 
half  throughout  these  results.  There  is  little  effect  on  the  particle  distribution  or 
potential  structure  if  the  temperature  ratio  is  changed  to  other  values.  The  electron 
temperature  is  expected  to  be  higher  than  the  ion  temperature  in  the  presheath  since 
electrons  absorb  energy  from  electric  and  magnetic  fields  faster  than  ions  and  other 
large  particles. 

Through  repeated  test  runs  of  the  model  it  was  found  that  fifty-one  points  in 
velocity  space  were  enough  to  provide  high  accuracy  and  produce  good  results. 
The  range  of  points  in  velocity  space  has  been  truncated  to  ±5  nondimens ional 
units.  The  results  show  that  at  ±5  the  distribution  is  near  zero,  substantiating  the 
truncation. 

A  time  step  of  0.2  nondimensionai  times  was  found  to  keep  the  solution  accurate. 
Three  nondimensionai  units  in  time  were  sufficient  to  produce  stable  results. 

It  was  found  that  the  magnitude  of  the  higher  order  terms  in  the  particle  velocity 
distribution  drop  off  very  rapidly.  Thus,  the  higher  order  terms  have  very  little 
impact  on  the  shape  of  the  potential  or  of  the  particle  distribution. 

To  understand  the  effects  of  a  quiescent  plasma  interacting  with  a  surface  the 
potential  gradient  in  the  neutral  plasma  has  been  set  to  zero.  To  accomplish  this 
the  a  term  in  the  expansion  of  potential  has  been  set  to  zero. 


U  =  Uq  +  aj  A  U  +  aiAU7  +  ••• 


(3.13) 
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where 

U0~az  and  At/  =  t**  (3.14) 

A  stable  solution  exists  only  for  a  specific  critical  value  of  the  exponential  coefi- 
cient,  p,  which  represents  the  scale  of  the  presheath.  The  quantitiy  p  is  nondimen- 
sionalized  as  <p0  -  >  where  n^r  is  an  ion  mean  free  path.  It  is  expected,  as  shown 

in  section  3.2,  that  the  critical  value  should  be  on  the  order  of  a  mean  free  path.  It 
was  found  that  *  0.4  produces  the  most  nearly  stable  results.  The  distributions 
become  unstable  for  values  greater  or  less  than  0.4.  The  small  remaining  instability 
at.  -A.  =  o.4  can  be  attributed  to  the  inexact  nature  of  the  numerical  solution. 

The  results  presented  here  are  first  order  and  produce  a  complete  picture  of  the 
structure  of  the  presheath  because  the  higher  order  terms  collective  contribution 
is  more  than  an  order  of  magnitude  smaller.  Figures  5.1  and  5.2  are  plots  of  the 
zeroth  and  first  order  expansions  of  the  particle  distribution  in  velocity  space.  The 
zero  order  term  remains  Maxwellian  because  the  potential  gradient  in  the  neutral 
plasma  is  zero.  The  first  order  term  of  the  distribution  obtains  a  profile  that  has 
roughly  the  shape  ( but  not  magnitude )  of  the  negative  first  derivative  of  the  zeroth 
order  solution.  The  potential  expansion  parameter  at  the  point  of  no  returning  ions 
is  determined  for  each  time  step.  Using  the  particle  distribution  functions  and  the 
known  potential  expansion  parameter  together  produce  the  overall  particle  velocity 
distribution  at  the  point  of  no  returning  ions,  the  presheath  -  sheath  interface. 
Figure  5.3  shows  this  distribution. 

The  positive  shift  in  the  total  distribution  is  as  expected  for  the  presheath.  The 
ions  are  being  pulled  into  the  wall.  The  particle  distribution  for  velocities  away  from 
the  wall  is  zero  for  the  case  of  no  returning  ions.  The  point  of  no  returning  ions 
exists  where  the  particle  distribution  for  velocities  away  from  the  wall  integrates 


7k 


33 


to  zero.  Figure  5.3  shows  that  the  ion  distribution  becomes  negative  for  velocities 
away  from  the  wall.  A  negative  particle  distribution  cannot  exist  physically.  The 
addition  of  higher  order  terms  does  not  correct  the  problem  because  the  expansion 
drops  off  so  quickly  that  any  higher  order  terms  have  no  impact  on  the  shape  of 
the  distribution.  The  problem  is  fundamental  to  the  type  of  collision  term  being 
applied  in  the  model.  The  Fokker  -  Planck  term  only  models  the  limit  of  small 
angle  collisions.  However,  large  angle  collisions  become  important  in  the  presheath. 

The  first  term  of  the  nondimensionalized  potential  structure,  <pi,  has  been  ini¬ 
tially  perturbed  to  1.0  x  io~4  to  obtain  the  results  presented  in  figures  5.1,  5.2  and 
5.3.  Perturbing  the  potential  structure  provides  the  model  with  the  nonequilibrium 
condition  necessary  to  initiate  the  time  dependent  development  of  the  presheath. 
The  strength  of  the  initial  perturbation  is  not  significant  to  obtaining  an  accurate 
particle  distribution  and  potential  structure  of  the  presheath.  Figures  5.4,  5.5,  and 
5.6  are  the  result  of  an  initial  perturbation  of  1.0  x  10~3  and  figures  5.7,  5.8  and  5.9 
are  the  result  of  an  initial  perturbation  of  1.0  x  10"5.  Comparing  these  results  show 
that  the  magnitude  of  the  initial  perturbation  only  affects  the  scale  of  the  first  order 
term  and  has  no  effect  on  the  overall  particle  distribution  in  velocity  space. 

Figure  5.10  is  a  plot  of  the  position  of  the  point  of  no  returning  ions,  the 
presheath  -  sheath  interface,  as  a  function  of  time  for  the  three  solutions.  Since  no 
source  of  ions  exists  in  the  model  the  relative  position  of  the  plasma  with  respect 
to  the  surface  changes  as  a  function  of  time.  The  wall  is  moving  into  the  plasma, 
or  the  plasma  is  moving  into  the  wall,  at  the  rate  at  which  the  wall  is  absorbing 
ions.  The  three  solutions  have  different  magnitudes  but  follow  the  same  profile. 
The  strength  of  the  perturbation  controls  the  relative  position  of  the  zero  point. 

Figure  5.11  is  a  plot  of  the  potential  structure  of  the  presheath  obtained  from 
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the  three  solutions  as  a  function  of  position.  The  data  for  the  potential  structure  has 
been  taken  from  the  solution  at  a  nondimensional  time  of  two.  Note  that  the  affect 
of  the  different  initial  perturbation  values  is  to  cause  a  shift  in  the  relative  position 
of  the  potential  but  has  no  effect  on  the  shape  of  the  potential  or  on  the  strength  of 
the  potential  at  the  point  of  no  returning  ions.  Changing  the  perturbation  strength 
alters  the  location  of  the  zero  point  but  not  its  shape.  The  stronger  the  perturbation 
the  further  the  zero  point  is  moved  from  the  surface.  The  horizontal  line  in  the  plot 
depicts  the  presheath  height  at  the  point  of  no  returning  ions.  The  vertical  lines 
show  the  position  of  the  point  of  no  .  ming  ions. 

A  tima  dependent  plot  of  the  presheath  height  at  the  point  of  no  returning  ions 
is  presented  in  figure  5.12.  This  plot  shows  that  the  time  evolution  of  the  sheath' 
height  approaches  smoothly  to  a  nearly  constant  value  of  0.16.  All  three  solutions 
fall  on  the  same  curve.  This  shows  that  the  strength  of  the  perturbation  does  no. 
affect  the  results  obtained. 
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CHAPTER  VI 
CONCLUSIONS 

The  solution  obtained  is  an  accurate  representation  of  the  time  dependent  de¬ 
velopment  of  the  Fokker  -  Planck  presheath.  The  model  produces  a  precise  potential 
structure,  however,  the  distribution  of  returning  ions  breaks  down  in  the  presheath. 
An  oscillation  developes  in  the  negative  tail  of  the  distribution,  as  seen  in  figures 
5.3,  5.6,  and  5.9.  This  oscillation  cannot  be  removed  by  including  additional  terms 
to  the  expansion.  In  addition,  the  sheath  height  of  0.16  determined  at  the  point  of 
no  returning  ions  is  roughly  an  order  of  magnitude  smaller  than  expected.  Both  of 
these  conditions  lead  to  the  conclusion  that  the  Fokker  -  Planck  collision  term  does 
not  represent  the  type  of  collisions  that  remove  the  returning  ions  in  the  presheath. 
This  breakdown  is  do  to  the  failure  of  the  Fokker  -  Planck  collision  term  to  model 
the  large  angle  collisions  that  take  place  within  the  presheath.  The  Fokker  -  Planck 
term  is  effective  at  modeling  the  collisions  present  in  the  center  of  the  plasma  but 
breaks  down  in  the  presheath.  The  primary  mechanism  behind  clearing  out  the 
returning  ions  from  within  the  presheath  is  not  particle  diffusion  as  represented  by 
small  angle  deflections  but  rather  the  large  velocity  changes  caused  by  large  angle 
collisions.  Since  the  Fokker  -  Planck  term  models  particle  collisions  that  represent 
the  limit  of  small  angle  collisions  it  is  inadequate  at  describing  the  mechanisms 
controlling  the  ion  velocity  distribution  moving  away  from  the  wall.  The  solutions 
obtained  using  a  Maxwellian  distribution  by  Bissell  and  Johnson!4)  and  Emmert  et 
al.t3l  and  those  obtained  using  a  charge  exchange  collision  model  by  RiemannM  and 
Main)7)  effectively  include  the  large  angle  collisions  since  they  model  the  collisions 
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by  instantaneous  changes  in  particle  velocity  and  position.  These  collision  models 
do  not  represent  the  Coulomb  collisions  taking  place  in  a  fully  ionized  plasma.  They 
do  not  represent  the  collision  processes  but  only  approximate  the  coilisional  effects. 

This  Fokker  -  Planck  presheath  model  produces  a  self-consistent  and  precise 
potential  structure.  The  particle  velocity  distribution  in  the  presheath  has  the 
correct  acceleration  of  ions  toward  the  wall  but  because  the  Fokker  -  Planck  collision 
term  only  models  the  limit  of  small  angle  collisions  it  is  unable  to  clear  the  particle 
distribution  of  returning  ions.  The  effect  of  not  modeling  the  large  angle  collisions 
b  that  the  particle  dbtribution  for  returning  ions  b  accurate  only  in  the  initial 
section  of  the  presheath  where  the  collbional  processes  are  dominated  by  particle 
diffusion.  The  collbional  processes  become  dominated  by  the  effects  of  the  large 
angle  collbions  as  the  interface  between  the  presheath  and  the  Debye  sheath  b 
approached.  Only  by  including  a  collision  term  which  accounts  for  these  large 
angle  collbions  can  a  presheath  model  produce  a  particle  velocity  dbtribution  that 
accurately  modeb  the  condition  of  no  returning  ions.  Thb  study  has  found  that  a 
presheath  model  which  describes  the  Coulomb  collbions  occurring  in  a  fully  ionized 
plasma  must  account  for  both  the  small  angle  collbions  and  the  large  angle  collisions. 
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Figure  3.1  Debye  Shielding 


Figure  3.2  Wall  Potential 


Presheath  Region(0(A,;)) 


Figure  3.3  Wall  Region  Model 
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Figure  4.1  Program  Diagram 


Figure  5.2  First  Order  Particle  Velocity  Distribution 
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Figure  5.6  Total  Particle  Velocity  Distribution,  pi  =  1.0  x  10' 
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Figure  5.7  Zero  Order  Particle  Velocity  Distribution,  =  1.0  x  10“ 


50 


i 


Nf(OIOtNON'f(00'-CS  +  (I!OCSN 
*-  OOOO  OOOO  I  CM 


till  I  I  I  I 

(9-301  a®uj.j_) 

(V)  tioi^nqinsiQ  aptixvj  «PJ0 


Figure  5.8  First  Order  Particle  Velocity  Distribution,  pi  =  10  x  10 
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Figure  5.10  Time  Dependent  Position  of  the  Potential 
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APPENDIX:  PROGRAM  LISTING 


THIS  PROGRAM  WAS  WRITTEN  BY  JEFFREY  P.  DANSEREAU 
THIS  VERSION  WAS  LAST  UPDATED  ON  7/17/87 


C 

C  THIS  PROGRAM  IS  A  TIME  DEPENDENT  MODEL  OF  THE  SHEATH  - 
C  PRESHEATH  OF  A  PLASMA.  IT  USES  A  FOKKER-PLANCK  COLLISION 
C  TERM  WITH  ONLY  COULOMB  COLLISIONS.  BELOW  IS  A  LIST  OF 
C  THE  VARIABLES  AND  THEIR  MEANING 
C 

C  T  -  DIAGONAL  REPRESENTATION  OF  T  MATRIX.  1ST  IS  DIAGONALS. 
C  2ND  IS  M  VEL  POS. 

C  V  -  3-D  MATRIX  OF  V  COMPONENTS.  1ST  POS.  IS  F’S,  2ND  IS 
C  N  PHI  POS.,  3RD  IS  M  VEL.  POS. 

C  W  -  3-D  MATRIX  OF  V  VALUES  AFTER  INVERSION  WITH  T  MATRIX. 
C  POS.  ARE  SAME  AS  V  MATRIX  WITH  ADDITIONAL  ROW  FOR  BC’S. 

C  D  -  2-D  MATRIX  OF  D  VALUES.  1ST  POS.  IS  F’S,  2ND  IS  M  VEL  POS. 

C  DD  -  2-D  MATRIX  OF  D  VALUES  AFTER  INVERSION  WITH  T  MATRIX 
C  PHI-  VECTOR  OF  PHI  VALUES 

C  F  -  2-D  MATRIX  OF  DENSITY  FUNCTIONS,  1ST  POS.  IS  THE  n  F’S 
C  THAT  ARE  BEING  USED.  (n-N-1).  F(1,X)«  FO  ECT...  THE 

C  2ND  POS  IS  THE  M  VEL  POS. 

C  TSTEP  -  VALUE  OF  DELTA  T  AS  TIME  IS  STEPPED  THROUGH 

C  VSTEP  •  VALUE  OF  DELTA  V  AS  VEL.  SPACE  IS  STEPPED  THROUGH 

C  RTIME  -  CURRENT  VALUE  OF  NON-DIM.  TIME 
C  ETA  -  CURRENT  VEL. 

C  NETA  -  ETA  PARAMETER  IN  POISSON  EQN 
C  TOTE  -  RATIO  OF  T  TO  T*(T  OVER  Tt) 


c 

SM 

-  SMALL  M(m)  IN  F-P  EQN 

c 

BM 

-  BIG  M(M)  IN  F-P  EQN 

c 

AA 

-  » IN  F-P  EQN 

C 

AP 

-  »’  IN  F-P  EQN 

C 

M 

-  NUMBER  OF  DIV.  IN  VEL.  SPACE 

c 

N 

-  NUMBER  OF  PHI  VALUES 

c 

ND 

-  NUMBER  OF  DIAGONALS  IN  T  MATRIX 

o  o 
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C  TIME  -  INTEGER  VALUE  IN  TIME  LOOP 
C  TEND  -  RTIME  TO  FINISH  SIMULATION 
C  L  -  INTEGER  TIME  VALUE  TO  END  SIMULATION 
C  FLAG1  -  FLAG  TO  PRINT  OR  NOT  PRINT  MATRICES (99  TO  PRINT) 

C  FLAG2  -  FLAG  TO  PRINT  OR  NOT  PRINT  INTERMEDIATE  MATRICES 

C  (99  TO  PRINT) 

C  TRBLE  -  FLAG  TO  PRINT  TROUBLE  STATEMENT  IF  MATRIX  IS  NOT 
C  INVERTABLE 

C  B.SOLN  -  INTERMEDIATE  VALUES  OF  VARIOUS  FUNCTIONS 
C  X,Y,Z  -  INTEGER  COUNTERS 
C 

C  THIS  PROGRAM  STEPS  THROUGH  TIME  SOLVING  THE  BOLTZMANN  EQN  FOR 
C  VALUES  OF  THE  EXPANDED  DENSITY  FUNCTION 
C 

REAL  T(6,8,202)IV(6,6,202),B(202),SOLN(202),TEND 
REAL  W(6,6,202)  ,D  (6,202)  ,DD(6, 202)  ,PHI(6)  ,L 

REAL  F(6,202),TSTEP,RTIME,SM,BM,AA,AP,VSTEP, ETA, NETA, TOTE 
REAL  FTOTAL(202),S0,DU1DUPHIl 

INTEGER  X,Y,Z,M)N,ND(TIME,TRBLE,FLAGl1FLAG2,FLAG3,PIPSTEPISKIP 

c 

READ  IN  PARAMETERS  AND  PRINT  THEM 

OPEN(UNIT-2,FILE=’INITDAT.DAT’,STATUS«’OLD’) 

READ(2,705)  MlNETA,TOTE,AA,AP,SM1FLAG3,SKIP,S0 
READ (2,707)  BM,VSTEP,N,ND,TEND,TSTEP,FLAG1,FLAG2 
CLOSE(UNIT*2) 

70S  F0RMAT(I4,1X,F6.4,1X,F6.4,1X,F6.4,1X,F6.4)1X,F6.4,1X,I3,1X, 

+  I3,1X,F9.4,/} 

707  F0RMAT(1X,F6.4,1X,F8.4,1X,I4,1X,I4,1X,F8.5I1X,F7.5,1X,I3,1X,I3) 
OPEN(UNIT-1,FILE=THLOUT\STATUS=,UNKNOWN’) 
OPEN(UNIT-3,FILE-’FPSHETH.OUT’1STATUS-,UNKNOWN’) 

WRITE(3,717)  . . 

+  '  . . 

717  FORMAT(A.A) 

WRITE(3,711)  ’FOKKER  -  PLANCK  SIMULATION  OUTPUT’ 

WRITE(3,717)  . . 

+  ,*•••••••••’ 

WRITE(3,710)  ’INPUT  PARAMETERS’ 

WRITE(3,715)  ’NUMBER  OF  VEL.  STEPS  -  M’,M 
WRITE(3,715)  ’NUMBER  OF  PHI  VALUES  -  N’,N 
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WRITE(3,715)  ’NUMBER  OF  DIAGONALS  IN  T  MATRIX  -  ND’.ND 
WRTrE(3,720)  ’VALUE  OF  ETA  IN  POISSON  EQN  -  NETA’.NETA 
WRITE(3,720)  ’VALUE  OF  T  OVER  TE  -  TOTE’, TOTE 
WRITE(3,720)  ’SIZE  OF  EACH  VEL.  STEP  -  VSTEP’.VSTEP 
WRITE(3,720)  ’SIZE  OF  EACH  TIME  STEP  -  TSTEP’.TSTEP 
WRITE(3,720)  ’VALUE  OF  ENDING  TIME  FOR  SIMULATION’, TEND 
WRITE(3,720)  *VALUE  OF  A  PARAMETER  IN  F-P  EQN  •  AA’.AA 
WRITE(3,720)  *VALUE  OF  A  PRIME  PARAMETER  IN  F-P  EQN  -  AP’.AP 
WRITE(3,720)  ’VALUE  OF  SMALL  M  IN  F-P  EQN  -  SM’,SM 
WRTTE(3,720)  ’VALUE  OF  BIG  M  IN  F-P  EQN  -  BM’.BM 
WRITE(3,716)  ’FLAG  TO  PRINT  PRIMARY  MATRICES(99  TO  PRINT)’, 

+  *  -  FLAGl’.FLAGl 

WRITE(3,718)  ’FLAG  MAKE  FDATA  AND  PHIDAT  NEW  FINAL  VALUES’, 

+  ’(99  *  YES)  -  FLAG3’,FLAG3 

WRITE(3,715)  ’TIME  SKIP  FOR  PRINT  OF  F  FILE  -  SKIP’, SKIP 
716  F0RMAT(SX,A,A,2X,I3) 

WRITE(3,710)  ’  ' 

710  FORMAT(15X,A,//) 

711  FORMAT(8X,A,/) 

715  FORMAT(5X,A,2X,I3) 

720  FORMAT(3X,A,2X,F9J) 

C 

C  MAKE  NON  TIME  DEPENDENT  QUANTITIES  AND  INITIAL  APPROXIMATIONS 
C  TO  F’S,  PHI’S,  AN  THE  BC 
C 

DU*0.0 
DUPHIl-0.0 
CALL  MAKEF(M.F) 

CALL  MAKEPHI(N,RTIME,PHI) 

CALL  TOT(F,FTOTAL,PHI,NIM,VSTEP,DU,DUPHIl) 

TT*0.0 

CALL  CONSERV(F,VSTEP,TT,M) 

C 

c  . . . 

c 

C  START  MAIN  TIME  LOOP 
C 


a  o 


58 


COMPUTE  END  LOOP  TIME  L 

L-TEND/TSTEP 
PSTEP-0 
C 

DO  40  TIME-1, INT(L)+1 

C  PRINT  RESULTS  OF  CURRENT  TIME  STEP 
C 

PSTEP— PSTEP+1 
IF  (TIME.LE.3)  GOTO  1000 
IF  (PSTEP.GE.SKIP)  THEN 
PSTEP— 0 
GOTO  1000 
ENDIF 
GOTO  1001 

1000  WRITE(3,350)  ’TIME  =  ’.RTIME.’TSTEP  -  ’.TIME-l 
DO  210  Z-l.N 

WRITE(3,250)  ’PHI(\Z.2,’)  -  ’,PHI(Z) 

210  CONTINUE 

WRITE(3,351)  ’DU  -  ’,DU,’DU*PHI(1)  -  ’.DUPHIl 
WRITE(3,711)  ’  ’ 

WRITE(3,275)  ’M’,’ETA’I,FO’1’Fl’,’F2’,’F3’,’FTOTAL’ 

ETA— 5.0 
DO  220  X-l.M 

WRITE(3,300)  X,ETA,F(l,X),F(2,X),F(3,X),FTOTAL(X) 

ETA-ETA+VSTEP 
220  CONTINUE 

1001  WRITE(3,710)  ’  ’ 

WRITE(1,813)  TIME, RTIME,PHI(1),PHI(2),?HI(3), DU, DUPHIl 
813  FORMAT(I3,1X,F7.4,1X,E13.6,1X,E13.0,1X,E13.0,1X,E13.0, 

+  1X.E13.8) 

IF  (TIME.GT.INT(L))  GOTO  40 

WRITER, 948)  ’CURRENTLY  IN  TIME  STEP’, TIME, ’OF’,INT(L) 

946  F0RMAT(5X,A,I4,2X,A,I4) 

C 

C  MAKE  TIME  DEPENDENT  QUANTITIES  T  MATRIX,  V  MATRIX,  D  MATRIX 
C 

DO  445  JJ— 1,2 

CALL  MAKET(ND,M,F,VSTEP,SM,BM,AA,AP,TSTEP,SO,PHI,N,T) 

CALL  MAKEV(N,M,F,VSTEP,V) 


I  01 


o  o 


CALL  MAKED(M,F,VSTEP)TSTEP1SM,BM,AA,AP1N>PHI,D) 
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INVERT  T  MATRIX  WITH  V  VECTORS,  MAKING  W  MATRIX 
C 

DO  50  X-l.N-l 
DO  60  Y=3,N 
DO  70  Z»1,M 

B(Z)s»V(X,Y,Z) 

70  CONTINUE 

CALL  MODIAG(M,ND,T,B,SOLN,TRBLE,X) 

IP  (TRBLE.EQ.999)  THEN 

WRITE(3,*)  ’MATRIX  HAS  NO  SOLUTION’ 
TRBLE=0 
END  IF 
DO  80  Z=1,M 

W(X,Y,Z)=SOLN(Z) 

80  CONTINUE 

60  CONTINUE 

50  CONTINUE 
C 

C  INVERT  T  MATRIX  WITH  D  VECTORS,  MAKING  DD  MATRIX 
C 

DO  90  X*1,N-1 
DO  100  Y=1,M 
B(Y)-D(X,Y) 

100  CONTINUE 

CALL  MODIAG(M,ND,T,B,SOLN,TRBLE,X) 

IF  (TRBLE.EQ.999)  THEN 

WRITE(3,*)  ’MATRIX  HAS  NO  SOLUTION’ 

TRBLE*0 
END  IF 

DO  110  Y=1,M 

DD(X,Y)-SOLN(Y) 

110  CONTINUE 

90  CONTINUE 

C 

C  GET  NEW  F  VALUES  FROM  NEW  PHI  VALUES 

C 

DO  160  X=1,N-1 

IF  ((JJ.EQ.1).AND.(X.EQ.2))  GOTO  160 


l  os 


IF  ((JJ.EQ.2).AND.(X.EQ.l))  GOTO  160 

DO  165  Z=1,M 

F(X,Z)=0.0 

-  1 

165 

CONTINUE 

DO  170  Y«3,N 

DO  180  Z=1,M 

F(X,Z)=sF(X,Z)+PHI(Y)*W(X,Y,Z) 

I 

180 

CONTINUE 

I 

170 

CONTINUE 

DO  190  Z=1,M 

I 

F(X,Z)=F(X,Z)+DD(X,Z) 

'  1 

190 

CONTINUE 

160 

C 

CONTINUE 

CALL  aVE{F,M,N,JJ) 

I 

445 

CONTINUE 

'  1 

C 

1 

c 

c 

c 

GET  NEW  PHI  VALUES 

CALL  DENSITY(N,M,PHI,F,TOTE,VSTEP} 

c 

c 

MAKE  TOTAL  DENSITY 

CALL  TOT(F1FTOTAL,PHI,N,M,VSTEPIDUfDUPHIl) 

I 

c 

c 

TT=RTIME+TSTEP 

CALL  CONSERV(F,VSTEP,TT,M) 

c 

c 

INCREASE  REAL  TIME  TO  NEXT  POSITION 

c 

RTIME-RTIME+TSTEP 

c 

c 

CONTINUE  TIME  LOOP 

40 

CONTINUE 

c 

c 

c 

c 

END  TIME  LOOP 

c 

* 

lOH- 

c 

C  FORMAT  STATEMENTS  FOR  PRINTS 
C 

250  F0RMAT(5X,A,I4)A,1X,E13.6) 

275  FORMAT{2XtA)2X,A19X,A,12X,A,12X,A,12X,A,10X,A) 

300  FORMAT(I3,lX,F5.2,lX,E13.6,lX,El3.6,lX,El3.8,lX  E13.6.1X.E12.5) 

350  FORMAT(5X,A,F7.4,10X,A,I5,/) 

351  F0RMAT(5X,A,E13.6,5X,A,E13.6) 

CLOSE(UNIT-3) 

IF  (FLAG3.EQ.99)  THEN 

OPEN(UNIT=4,FILE«’FDATA.DAT\STATUS-,UNKNOWN’) 
DO  265  X-1,M 

WRITE(4,266)  F(1,X),F(2,X)1F(3,X)IF(41X) 

266  FORMAT(El3.8,lX,El3.6,lX,El3.6,lX,El3.6) 

265  CONTINUE 

close(unit-i) 

OPEN(UNIT»8,FILE»,PHIDAT.DAT,,STATUS»,UNKNOWN') 
DO  267  X-1,N 

WRITE(8,266)  RTIME 
WRrrE(8,268)  PHI(X) 

268  F0RMAT(E13.6) 

267  CONTINUE 
END  IF 

STOP 

END 


REAL  F(6,202),ETA,VSTEP,NETA)TOTE,B,PI.C 
INTEGER  X,M,Z 

OPENfUNrr^'S.Fn.EaB’FDATA.DAT’.STATUS^’UNKNOWN’) 
ETA— 5.0 

PRINT*, ’INPUT  M,VSTEP’ 

READ(*,*)  M.VSTEP 
N«M-(M-l)/2 
DO  10  X-1.N 

F(1,X)-EXP(-(ETA**2)) 

F(2,X)»0.0 

F(3,X)»0.0 

F(4,X)=*0.0 

ETA-ETA+VSTEP 


10  CONTINUE 
Z-l 

DO  30  X=M,N+l,-i 
DO  40  Y=l,4 

F(1.X)»F(1,Z) 

40  CONTINUE 

Z-Z+l 

30  CONTINUE 
DO  20  X*1,M 

WRITE(3,100)  F(l,X),F(2,X),F(3,X),F(4fX) 
100  FORMAT(El3.6,  lX.El3.fi,  lX.El3.fi,  IX, El3.fi) 

20  CONTINUE 
STOP 
END 


REAL  PHI(6),RTIME 
INTEGER  N.X 
N=5 

RTIME— 0.0 

WRITE(V)  -INPUT  PHI(-l)’ 

READ(*,200)  PHI(l) 

WRITEC, *)  ’INPUT  PHI(O)' 

READ(*,200)  PHI(2) 

200  FORMAT(F10.5) 

WRITEC,*)  ’INPUT  PHI(l)’ 

READ(*,200)  PHI(3) 

DO  10  X-4,N 
PHI(X)-0.0 
10  CONTINUE 

OPEN(UNIT*9,FILE=’PHIDAT.DAT’, STATUS— ’UNKNOWN’) 
WRITE(9,120)  RTIME 
DO  20  X-1,N 

WRITE(9,120)  PHI(X) 

120  F0RMAT(E13.6) 

20  CONTINUE 
RETURN 
END 
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SUBROUTINE  CONSERV(F,VSTEP,T,M) 

REAL  F(8,202),F0(202),  VSTEP, ZF(202),Z2F(202),ENER, MOM, DEN, T, ETA 
INTEGER  X,M 
ETA=-5.0 
DO  10  X=1,M 
F0(X)-F(1,X) 

ZF(X)-F(1,X)*ETA 
Z2F(X)*ZF(X)*ETA 
ETA-ETA+VSTEP 
10  CONTINUE 

CALL  SIMPS(F0,M, VSTEP, DEN) 

CALL  SIMPS(ZF,M,VSTEP,MOM) 

CALL  SIMPS(Z2F,M, VSTEP, ENER) 

OPENtUNIT-O.FtLEa’CONSERV.OUT’.STATUS^'UNKNOWN’) 
WRITE(9,100)  T,DEN,MOM,ENER 
100  FORMAT(2X,F8.S,lX,El3.8,lX,El3.8,lX,El3.8) 

RETURN 

END 


SUBROUTINE  DENSITY(N,M, PHI, F, TOTE, VSTEP) 

REAL  F(6, 202), PHI(8), TOTE, VSTEP, G(202),R(4),NO(6) 

INTEGER  M.N.X.Y 
DO  10  X-1,N-1 
DO  20  Y=l,M 
G(Y)-F(X,Y) 

20  CONTINUE 

CALL  SIMPS(G,M,  VSTEP, NO(X)) 

10  CONTINUE 

R(l)«NO(2)/NO(l) 

R(2)-NO{3)/NO(l) 

R(3)«NO(4)/NO(l) 

PHI(3)— PHI(2)*R(1)/(2.0*TOTE) 

Pffl(4)-PHI{2)*((((PHI(3)/PHI(2))**2)*((2.0*TOTE)**2)/2.0)-R(2)) 
+  /(2.0'TOTE) 

PHI(5)»PHI(2)*((((2.0*TOTE)**2)*PHI(3)*PHI(4)/(PHI(2)**2)) 

+  -(((2.0*TOTE)**3)*((PHI(3)/PHI(2))**3)/6.0)-R(3))/(2.0*TOTE) 

RETURN 

END 
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SUBROUTINE  MAKEF(M.F) 

REAL  F(6,202) 

INTEGER  X,M 

OPEN(UNIT=fl,FILw=,FDATA.DAT’,STATUS=’UNKNOWN’) 
DO  10  X=1,M 

READ(S.IOO)  F(1,X),F(2,X),F(3,X),F(4,X) 

100  F0RMAT(E13.6,1X,E13.6,1X,E13.6,1X,E13.6) 

10  CONTINUE 
RETURN 
END 


SUBROUTINE  MAKEPHI(N,RTIME,PHI) 

REAL  PHI(6),RTIME 
INTEGER  N,X 

OPEN(UNIT*9IFILE=,PHIDAT.DAT^STATUS=■,UNKNOWN,) 
READ(9,120)  RTIME 
DO  20  X*1,N 

READ  (9, 120}  PHI(X) 

120  F0RMAT(E13.6) 

20  CONTINUE 
RETURN 
END 


SUBROUTINE  MAKEV(N,M,F,VSTEP,V) 

REAL  F(fl, 202), V(6, 6,202), VSTEP,FD1(6),Z, ETA 
INTEGER  N,M,X,Y,R,S 
ETA— 5.0 
DO  10  X-1,M 

CALL  FDl(F,N,M,X,VSTEP,FDl) 

IF  ((X.EQ.l).OR.(X.EQ.M))  THEN 
DO  15  Y-l.N-l 
DO  17  R*3,N 
V(Y,R,X)=0.0 
17  CONTINUE 

15  CONTINUE 

ELSE 

DO  20  Y*1,N-1 
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Z=1.0 

S=1 

DO  30  R=3,N 

IF  (S.GT.(Y-l))  THEN 
V(Y,R.X)=0.0 
ELSE 

V(Y,R,X)=*-Z*FDl(Y-R+2) 
END  IF 
S-S+l 
Z=Z+1.0 

30  CONTINUE 

20  CONTINUE 

ENDIF 

ETA=ETA+VSTEP 
10  CONTINUE 
RETURN 
END 


SUBROUTINE  MAKEDfM.F.VSTEP.TSTEP.SM.BM.AA.AP.N.PHI.D) 
REAL  VSTEP.TSTEP.SM, BM,AA,AP,F(6, 202), D(6, 202)^(6,202), PHI(6) 
REAL  FD1(6)IFD2(6),AH(202),BH(202),B(6,202) 

REAL  ETA,GDl,GD2 
INTEGER  N,M,X,Y,Z 
DO  10  X-l.N-i 

CALL  GETAB(F,VSTEP,M,SM,BM,AA,AP,AH,BH,X) 

DO  15  Y=1,M 
A(X,Y)=AH(Y) 

B(X,Y)=BH(Y) 

15  CONTINUE 
10  CONTINUE 
DO  17  Y-l.N-l 
D(Y,M)«0.0 
D(Y,1)»0.0 
17  CONTINUE 

ETa— 5.C+VSTEP 
DO  20  X«2  M-1 
DO  V  •  -1.N-1 
7  X)-0.0 
L’O  Z=*1,Y 
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IF  (Z.LE.Y-1)  THEN 

GD2=(A(Z+l,X+l)*F(Y-Z,X+l)-2.0* 

+  A(Z+1,X)*F(Y-Z,X)+A(Z+1,X-1)*F(Y-Z,X-1))/(VSTEP**2) 

GDl-(B(Z+l,X+l)*F(Y-Z,X+l)-B(Z+l 
+  ,X-1)*F(Y-Z,X-1))/(2.0*VSTEP) 

D(Y,X)=D(Y,X)+GD2+GDl 
END  IF 

40  CONTINUE 

D(Y,X)*D(YIX)+F(Y1X)/TSTEP 
30  CONTINUE 

ETA-ETA+VSTEP 
20  CONTINUE 
RETURN 
END 


SUBROUTINE  FDl(F,N,M,X, VSTEP, FDl) 

REAL  F (6,202)  ,FD  1  (6) , VSTEP 
INTEGER  Y,X,N,M 
DO  10  Y-l.N-1 

IF  (X.LE.2)  THEN 

FDl(Y)=(-F(Y1X+2)+4.0*F(Y,X+l)-3.0*F(Y,X))/(2.0*VSTEP) 
ELSE  IF  (X.GE.M-1)  THEN 

FDl(Y)=*(3.0*F(Y,X)-4.0*F(Y,X-l)+F(Y,X-2))/(2.0*VSTEP) 

ELSE 

FDl(Y)-{F(YIX+l).F(Y,X.l))/{2.0' VSTEP) 

ENDIF 

10  CONTINUE 
RETURN 
END 


SUBROUTINE  MAKET(ND,M,F,VSTEPISM,BM,AA,AP,TSTEP,SO,PHI1N,T) 
REAL  VSTEP,TSTEPlSM,BM,AA,AP1F(e,202),T(6,8,202),A(202) 

REAL  ETA,B(202),PHI(6),P 

INTEGER  ND.M.X.Y.N.Z 

CALL  GETAB(F, VSTEP, M)SM,BM,AA,AP,A,B,1) 

P-0.0 

DO  5  Z-l.N-l 
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ETA— 5.0 
DO  10  X«1,M 

IF  ((X-EQ.M).OR.(X.EQ.l))  THEN 
T(Z,2,X)=0.0 
T(Z,3,X)=1.0 
T(Z,4,X)=0.0 
ELSE 

T(Z,2,X)=(((-A(X-l)/VSTEP)+B(X-l)/2.0)/VSTEP) 

T(Z,2,X)*T(Z,2,X)-PHI(1)/(2.0*VSTEP) 

T(Z)3,X)»1.0/TSTEP+(2.0*A(X)/(VSTEP**2)) 

T(Z,3,X)=T(Z)3,X)+PHI(2)'P*ETA 

T(Z,4,X)=((-A(X+l)/VSTEP-B(X+l)/2.0)/VSTEP) 

T(Z,4,X)-T(Z,4,X)+PHI(1)/(2.0*VSTEP) 

ENDIF 

ETA-ETA+VSTEP 
10  CONTINUE 

P=P+1.0 
5  CONTINUE 
RETURN 
END 


SUBROUTINE  GETAB(F,VSTEP,M,SM,BM,AA>AP,A,B,Y) 
REAL  A(202)IB(202)lAA,AP,F(6,202),VSTEP,Z,SMtBM,E(202) 
REAL  P(202),H(202),G(202),K(202),DB,DA,Sl,S2,S3,S4,S5 
INTEGER  MAY.N 
Z— 5.0 

DO  10  X-1,M 

CALL  FIND A(F VSTEP, AA, AP,M,A(X),Y) 

CALL  FINDB(F,Z1VSTEPIAA,AP,SM,BM,M,B(X)fY) 
Z-Z+VSTEP 
10  CONTINUE 
ETA— 3.0 
DO  20  X-1,M 

H(X)-B(X)*F(Y,X) 

G(X)-A(X)*F(Y,X) 

K(X)-ETA*B(X)*F(Y,X) 

P(X)=ETA*F(Y,X) 

E(X)-F(Y,X) 

ETA-ETA+VSTEP 
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20 


30 


CONTINUE 

CALL  SIMPStH.M.VSTEP.Sl) 
CALL  SIMPS(G,M,VSTEP,S2) 
CALL  SIMPS(K,M,VSTEP,S3) 
CALL  SIMPS(P1M,VSTEP,S4) 
CALL  SIMPS(E,M,VSTEPISS) 
IF  (SS.EQ.0.0)  RETURN 
DB— S1/S5 

DA— (S2-S3-DB'S4)/S5 
DO  30  X-1,M 

B(X)=B(X)+DB 

A(X)*A(X)+DA 

CONTINUE 

RETURN 

END 


SUBROUTINE  FIND A(F ,Z,VSTEP,AA,AP,M, A,Y) 

REAL  F(6,202),2,VSTEP,AA,AP,A,ETA,SOLN,H(202),FD2 
INTEGER  X,M,Y 
ETA— 5.0 
DO  10  X-1.M 

CALL  FD2(FIM,X,VSTEP,FD2,Y) 
H(X)*FD2*ABS(ETA-Z) 

ETA-ETA+VSTEP 
10  CONTINUE 

CALL  SIMPS(H,MlVSTEP,SOLN) 

A»(SOLN*AA)/(2.0*AP) 

RETURN 

END 


SUBROUTINE  FINDBJF.Z.VSTEP.AA.AP.SM^M^^Y) 

REAL  F(6,202)IZ,VSTEP,AA,APIB,ETAIG2JSOLN 

REAL  J(202),FD3iSM)BM 

INTEGER  X,M,Y 

ETA— 5.0 

DO  10  X-1.M 

CALL  FD3(F)M,X,VSTEP1FD3,Y) 
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J(X)=FD3*ABS(ETA-Z) 
ETA=ETA+VSTEP 
10  CONTINUE 

'CALL  SIMPS(J,M,VSTEP,SOLN) 
B— (AA/AP)*(SM/(2.0*BM))*SOLN 
RETURN 
END 


SUBROUTINE  SIMPS(F,N,H, RESULT) 

REAL  F(202),H, RESULT 

INTEGER  N.NPANEL.NHALF.NBEGIN.NEND 

NPANEL=*N-l 

NHALF=NPANEL/2 

NBEGIN-i 

RESULT=0.0 

IF  ((NPANEL-2*NHALF).NE.O)  THEN 

RESULT=»3.0*H/8.0*(F(l)+3.0*F(2)+3.0*F(3)+F(4)) 
NBEGIN-4 
IF  (N.EQ.4)  RETURN 
ENDIF 

RESULT=RESULT+H/3.0*(F(NBEGIN)+4.0*F(NBEGIN+1)+F(N)) 

NBEGIN-NBEGIN+2 

IF  (NBEGIN.EQ.4)  RETURN 

NEND-N-2 

DO  10  I«NBEGIN,NEND,2 

RESULT«RESULT+H/3.0*(2.0'F(I)+4.0*F(I+1)) 

10  CONTINUE 
RETURN 
END 


SUBROUTINE  FD2(F,M,X,VSTEP,FD2,Y) 

REAL  F(«,202),FD2,VSTEP 
INTEGER  X,M,Y 
IF  (X.LE.1)  THEN 

FD2*(F(Y,X+2)-2.0*F(Y,X+l)+F(Y,X))/(VSTEP**2) 
ELSE  IF  (X.GE.M)  THEN 

FD2»(F(Y,X).2.0*F(YlX-l)+F(Y,X-2))/(VSTEP**2) 
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ELSE 

FD2=(F(Y,X+1)-2.0*F(Y,X)+F(Y,X-1))/(VSTEP*,2) 
END  IF 
RETURN 
END 


SUBROUTINE  FD3(F,M,X,VSTEP,FD3,Y) 

REAL  F(8,202),FD3,VSTEP 
INTEGER  X,M,Y 
IF  (X.LE.2)  THEN 

FD3*(F(Y^C+3)-3.0*F(Y,X+2)+3.0*F(Y,X+l)-F(Y,X))/(VSTEP**3) 
ELSE  IF  (X.GE.M-1)  THEN 

FD3*s(F(Y,X)-3.0*F(Y,X-l)+3.0*F(Y,X-2)-F(Y,X»3))/(VSTEP**3) 

ELSE 

FD3»(F(Y1X+2).2.0*F(Y1X+l)+2.0*F(Y,X-l)-F(Y1X-2))/ 

+  (2.0*(VSTEP**3)) 

ENDIF 

RETURN 

END 


SUBROUTINE  MODIAG(M,D,MATRIX,C,SOLN1TRBLE,N) 
REAL  A(8,202),B,SOLN(202),MATRIX(618,202),C(202) 
INTEGER  M,D,W,X,Y,Z>TRBLEIN,IIMNl 
TRBLE-0 
DO  10  X=l,D+2 
DO  20  Y=1,M 

IF  (X.EQ.D+2)  THEN 
A(X,Y)-C(Y) 

ELSE 

A(X,Y)=*MATRIX(N)XtY) 

ENDIF 

20  CONTINUE 
10  CONTINUE 
DO  30  I-2,M 

A(2,I)-A(2,I)/A(3,I-1) 

A(3,I)-A(3,I)-A(2,I)*A(4II-1) 

A(5,I)-A(5,I)-A(2,I)*A(5,M) 
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30  CONTINUE 
DO  80  X=1,M 

IF  (A(3,X).EQ.O)  THEN 

PRINT*, ’MATRIX  HAS  NO  SOLUTION’ 
TRBLE=999 
GOTO  999 
END  IF 

80  CONTINUE 
MNl-M-1 

A(5,M)=sA(5,M)/A(3,M) 

DO  40  I»MN1,1,-1 

A(5,I)-(A(5,I)-A(4,I)*A(5,I+1))/A(3,I) 

40  CONTINUE 
DO  50  X-l.M 

SOLN(X)=A(5,X) 

50  CONTINUE 
999  RETURN 
END 


SUBROUTINE  AVE(F,M,N,JJ) 

REAL  F(6,202) 

INTEGER  X,N,M,Z,JJ 
DO  443  Z-l.N-l 

IF  ((JJ.EQ.1).AND.{Z.EQ.2))  GOTO  443 
IF  ((JJ.EQ.2).AND.(Z.EQ.l))  GOTO  443 
DO  444  X-2.M-1 

F(Z,XH.025*F(Z)X-1)+F(Z,X)+.025*F(Z,X+1))/1.05 
444  CONTINUE 
443  CONTINUE 
RETURN 
END 


SUBROUTINE  TOTfF.FTOTAL.PHLN.M.VSTEP.DU.DUPHIl) 
REAL  F(8,202)1FTOTAL(202),PHI(6),VSTEP,G(202)1NO(«)IDU 
REAL  Xl,X2,TOL,Fl,F2,XERR 
INTEGER  M,N,X,Y,NLIM 
L-(M+l)/2 
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DO  10  X-l.N-l 
DO  20  Y«1,L 
G(Y)-F(X,Y) 

20  CONTINUE 

CALL  SIMPS(G,L,VSTEP,NO(X)) 

10  CONTINUE 
Xl-0.0 
X2-1.0E0 
TOL-.0001 
NLIM-50 

30  CALL  FCN(NO,XltN,Fl) 

CALL  FCN(NO,X2,N,F2) 

IF  (Fl*F2.GT.0.0)  THEN 
X2=X2*10.0 

IF  (X2.GT.1.0E1S)  THEN 

WRITE(*,*)  ’NO  SIGN  CHANGE  UP  TO  1E15’ 
RETURN 
END  IF 
GOTO  30 
ENDIF 

DO  40  J»1,NLLM 
DU-(Xl+X2)/2.0 
CALL  FCN(NO,DUtN,FR) 

XERR=*ABS(Xl-X2)/2.0 
IF  (XERR.LE.  TOL)  GOTO  1000 
IF  (ABS(FR).LE-TOL)  GOTO  1000 
IF  (FR‘XI.GT.0.0)  THEN 
X1»DU 
F1«FR 
ELSE 

X2-DU 

F2-FR 

ENDIF 

40  CONTINUE 

WRITE(V)  ’NLIM  EXCEEDED’ 

RETURN 

1000  DOS0X-1.M 

FTOTAL(X)«0.0 
DO  60  Y-l.N-l 

FTOTAL  (X)  »FTOTAL  (X) +DU*  *  (Y-l)  *F(Y,X) 
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APPENDIX  B  -  The  TEC  program 
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noon 


C  DOS  FILE  TECINIT .FOR 

C* ***************************************************************** r 

TEC  DATA  INITIALIZATION  ACCESS  CODE 
WRITTEN  BY  GREGORY  L  RIDDERBUSCH 
AUGUST  1986 

****************************************************************** 
REAL  ARRAY (9) 

INTEGER  IPARAM (5) 

WRITE (*,103) 

OPEN (UNIT-2,  FILE-' PINDATA.DAT/ , STATUS-' OLD' ) 

READ (2/101)  (ARRAY (I) / 1=1/  9) 

READ (2/102)  (IPARAM (I) ,1=1,5) 

10  WRITE (*,104)  (ARRAY (I) ,1-1,7) 

WRITE (*,  105) 

READ (*, *)  I VALUE 
IF  (I VALUE  .EQ.  0)  THEN 
GOTO  20 
ELSE 

WRITE (*,106) 

READ (*, *)  VALUE 
ARRAY  ( I VALUE ) -VALUE 
GOTO  10 
END  IF 

20  WRITE (*,107)  (IPARAM (I) ,1=1,5) 

WRITE (*,110)  ARRAY (8) , ARRAY (9) 

WRITE (*,105) 

READ (*, *)  IVALUE 
IF  (IVALUE  .EQ.  0)  THEN 
GOTO  30 
ELSE 

WRITE (*, 108) 

IF  (IVALUE  .EQ.  6  .OR.  IVALUE  .EQ.  7)  THEN 
READ (*, *)  VALUE 
ARRAY ( IVALUE +2 ) -VALUE 
GOTO  20 

FT  QF 

READ  (*, *)  I NEW 
IPARAM (IVALUE) -INEW 
GOTO  20 
END  IF 
END  IF 

30  REWIND (2) 

WRITE (2,101)  (ARRAY (I) ,1-1, 9) 

WRITE (2,102)  (IPARAM(I) , 1-1,5) 

CLOSE (2) 

WRITE (*,109) 

STOP 

101  FORMAT (F8 .  1/F8 . 1/F6 . 3/F6 . 3/F6 . 3/F6 . 3/F7 . 2/F5 .2/F6 . 1) 

102  FORMAT (I1/I1/I3/I3/I3) 

103  FORMAT (IX, '  ******************************************/ / 

S'  TEC  OPERATING  CONDITIONS'/ 

SIX, '******************************************>) 

104  FORMAT (//4X, 'CURRENT  CONVERTOR  OPERATING  SETTINGS:'// 


S7X, '1. 

EMITTER  TEMPERATURE: 

' , F8 . 1, '  KELVIN',/ 

S7X, '2. 

COLLECTOR  TEMPERATURE: 

' , F8 . 1, '  KELVIN',/ 

S7X, ' 3 . 

EMITTER  WORK  FUNCTION: 

' , F6 . 2, ' 

EV'  / 

S7X, '4. 

COLLECTOR  WORK  FUNCTION: 

' ,  F6 . 2, ' 

EV'/ 

S7X, '5. 

CONVERTOR  PRESSURE: 

' ,F6.2, ' 

TORR' / 

S7X, '6. 

GAP  THICKNESS: 

' ,F6.2, ' 

MM'/ 

il* 


47X, '  7 «  OPERATING  CURRENT:  ',F7.2,'  AMPS/CMA2') 

105  FORMAT (/AX, 'ENTER  ID  NUMBER  OF  VALUE  TO  BE  CHANGED,  0=NONE:  ') 

106  FORMAT (4X, 'ENTER  NEW  OPERATING  SETTING:  ') 

107  FORMAT (//4X, 'CURRENT  TEC  FUNCTION  SETTINGS:'// 

47X, '1.  DIAGNOSTIC  LEVEL:  ',13,'  (0-RESTRICTED  OUTPUT)'/ 

47X, '  ( 1 -FULL  OUTPUT) ' / 

47X, '  (2-ENABLE  SHEATH  DIAGNOSTICS)'/ 

47X, '  (3-ENABLE  DOT  DIAGNOSTICS)'/ 

47X, '2.  RESTART  SEQUENCE:  ',13,'  (0-DEFAULT  STARTUP  VALUES)'/ 
47X, '  (1-RESTART  WITH  PREVIOUS  VALUES)'/ 

47X, '3.  STEPS  BETWEEN  PRINTS:  ',13/ 

47X,'4.  POINT  DENSITY:  ',13,'  (11, 21, 31,  . . . 151) ' / 

47X, ' 5 .  LOTUS  SKIP  FACTOR:  ',13,'  (1..99)') 

108  FORMAT (4X, 'ENTER  NEW  FUNCTIONAL  SETTING:  ') 

109  FORMAT (IX, '******************************************' / 

4'  '/ 

41X, '******************************************') 

110  FORMAT (7X, '6.  TIME  STEP:  ',F5.2, 

4'  (0.1, 0.2, 0.3, 0.4, 0.5) ' /7X, ' 7.  END  TIME:  ', 

4F6.1,'  (1.0, 2.0, ...,10.0)') 

END 


C  DOS  FILE  PRED1 . FOR 

C**********************************************************************!H r** 

C  PROGRAM  TEC 

q*  ***  **** *********  ****** *  **********  ************  **********  *  *  *  ********  *  *  *  *  *  *  *  * 
REAL  CNE, ENE, ECHI, CCHI, EALP HA, CALPHA, LAMNEB, LAMTAU 
REAL  DTP, T2, AN/ AT, CN, CT, BN, BT, RE, KN/ TCHAR, PN, DT 
REAL  SMR/ LAMDAR/NR/ TE, TC, ENR, I , ARECN, EGNDB, ELOSSB, NNR 
REAL  TAD (0 : 150) /NEB (0 : 150) / DELTAT/ SN, ST, PI, CA, CSAHA, DZ 
REAL  DTAUNDZ/MUI (0:150) , RMUR, TAUN (0 : 150) , EMISS, TIME2, LCCHI 
REAL  NDOT1 (0 : 150)  /  TDOT1 (0:150) , NNB (0 : 150) / TIME1, JNET/ CV  (0 : 150) 

REAL  FYEN/ IVD/ IKN, EGRADA/ PHIB, EWF, CWF/ D/ESOURCE  (0:150) 

INTEGER  N,  IDEN, EFIX, CFIX, CHKDOT, ICOUNT/ NSTEPS/ C, PC/ LS/ LC, EO, EC 
C 

COMMON  /PRED/  CA/ CSAHA/DT/DTAUNDZ/DZ/ EGNDB, ELOSSB, ENR/ I, IDEN/KN 
COMMON  /PRED/  LAMDAR, LAMNEB, LAMTAU, MUI , N,NNB, NNR, NR, PI, RE, RMUR 
COMMON  /PRED/  TAUN, EFIX, CFIX, CHKDOT, FYEN, ENE, ECHI, CCHI, CNE 
COMMON  /ERED/  EALP HA, CALPHA, CV, E SOURCE, LCCHI 
C 

data  nsteps  /I / 

C 

OPEN  (UNIT-2,  FILE-' PINDATA.DAT' , STATUS*' OLD' ) 

OPEN  (UNIT-4, FILE*' EXOUT.DAT' , STATUS*" UNKNOWN' ) 

OPEN  (UNIT-7,  FILE*' LOTUSl.DAT' /STATUS-' UNKNOWN' ) 

OPEN  (UNIT-8, FILE*' PREDRES . DAT' , STATUS*' OLD' ) 

OPEN  (UNIT-9, FILE*' PRNTOUT.DAT' /STATUS-' UNKNOWN' ) 

C _ SET  ABBREVIATED  PRNTDAT  SAMPLING  POINTS  WITH  PC 

C...SET  LOTUS  SKIP  COUNTER  WITH  LC 
C 

CALL  INITIAL (TE,TC, EWF, CWF, PN, NSTEPS, DTP, T2, AN, AT, 

&CN, CT,BN, BT, TCHAR, SMR, ARECN, DELTAT, SN, ST, TAU, NEB, LS, pc) 
c-pc 
LC-LS 
EO-100*LS 
EC-EO 
C 

IF  (CHKDOT  .EQ.  3) 

&  OPEN  (UNIT-9, FILE*' D:PDOTDIAG. DAT' , STATUS-' UNKNOWN' ) 

IF  (CHKDOT  .GT.  1) 

&  OPEN  (UNIT-11, FILE«'D:PSTHDIAG. DAT', STATUS-' UNKNOWN' ) 

C 

C - 

EGRADA- 0 . 0 

ICALCS-INT (T2 /DTP+0 .001) 

WRITE (*,76) 

DO  30  ICOUNT-0, I CALCS 

WRITE  (*,77)  (icount) ,100* (icount) /ICALCS 
TIMEl-(icount) *DTP 
TIME2- (icount+1) *DTP 

CALL  PREDCOR(TIMEl, TIME2, TAU, NEB, NSTEPS, TE, TC, TDOTl,NDOTl, 

6  PN,  ARECN,  TCHAR,  AN,  AT,  BN,  BT,  CN,  CT,  IVD,  SMR,  EGRADA,  PHIB,  JNET) 

VOUT-EWF  +  (PHIB  +  IVD /I  +  LCCHI) *TE/11600  -  CWF 
IF  (C  .EQ.  PC)  THEN 

WRITE  (9,67)  TIME1, VOUT 
IF  (CHKDOT  .EQ.  0)  THEN 
WRITE  (9,70) 

WRITE  (9,71)  (II, NEB (II) , TAU (II) , 11*1, N) 

ELSE 

EMISS-ENR/NEB (1) 

WRITE (9,  72)  ENE, ECHI, EALP HA, CNE, LCCHI +CCHI , CALPHA, PHIB, 


£  IVD/I/ EMISS 

WRITE (9,73)  (I1,ND0T1 (II) ,  NEB (II)  , 

£  TDOT1 (II)  , TAU (II)  , 11=0, N+l) 

END  IF 

IF  (EFIX  .EQ.  1)  WRITE (9, 74) 

IF  (CFIX  .EQ.  1)  WRITE (9/ 75) 

00 
END  IF 
C-C+l 

IF  (LC  .EQ.  LS)  THEN 

WRITE (7, 78) TIME1, VOUT, ENE, JNET, ECHI, CCHI, EALPHA, CALPHA, 

£  PHIB, IVD/I, EMISS, NEB (1) , NEB (11) , TAU (1) , TAU (11) 

LOO 
END  IF 
LC-LC+1 

IF  (EC  .EQ.  EO)  THEN 
WRITE (4,*)  TIME1 

WRITE (4, 71)  (II, NEB (II) , TAU (II) ,11-1, N) 

EOO 
END  IF 
EOEC+1 

30  CONTINUE 

C - - - 

c*****0(JTPUT  STARTUP  VALUES  TO  PREDRES.DAT 

WRITE (8,  69)  ENE, CNE, ECHI, CCHI, EALP HA, CALPHA,N 
WRITE (8,68)  (NEB (II) , TAU (II) , I 1-0, N+l) 

REWIND (8) 

CLOSE  (8) 

STOP 

C 

67  FORMAT  ( /  / , ' RESULTS  AT  TIME  -  ' , F8.2, 4X, 'OPERATING  VOLTAGE-' , F6 . 3) 

68  FORMAT (2F8. 3) 

69  FORMAT (F9.4/F9.4/F8.3/F8.3/F9.4/F9. 4/13) 

70  FORMAT (/3X,'  #  NEB (#)  TAU(#)'/ 

£3X, ' - ') 

71  FORMAT (3X,I3,4X,F6.3,7X,F5.2) 

72  FORMAT (/3X, 'ENE  -' , F8 . 3, 4X, ' ECHI-' , F8 . 3, 4X, 'EALPHA-' , F8 . 3/ 

&3X, ' CNE  -' ,F8 . 3, 4X, ' CCHI-' , F8 . 3, 4X, ' CALPHA-' , F8 . 3/ 
£3X,'PHIB-',F8.3,4X,'VD  -' , F8 . 3, 4X, ' EMISS  -',E8.3/ 


£/3X, '  I  NDOT (#)  NEB (#)  TDOT (#)  TAU(#)', 

£/3X,  ' - ' ) 


73  FORMAT (3X, 13, 4X, F7 . 4,  4X, F6 . 3, 7X, F7 . 4, 4X, F5 . 2) 

74  FORMAT ( ' *  *  *AT  LEAST  ONE  UNPHYSICAL  EMITTER  BC  WAS  INVOKED.') 

75  FORMAT C'***AT  LEAST  ONE  UNPHYSICAL  COLLECTOR  BC  WAS  INVOKED.') 

76  FORMAT (/IX,  '  TEC  CALCULATIONS  BEGIN. . . (WAIT) . . . ' ) 

77  FORMAT  (4X,' ITERATION  #',I3,'  COMPLETION— ',  13, '  %' ) 

78  FORMAT (16E13. 6) 

END 

£*************************************************************************** 
SUBROUTINE  PREDCOR (T1 , T2 , TAU, NEB, NSTEPS, TE, TC, TDOT 1, NDOT 1, 

+  PN, ARECN, TCHAR, AN, AT, BN, BT, CN, CT, IVD, SMR, EGRADA, PHIB, JNET) 

£*  *************************  ■  v  ************************************************ 

REAL  Tl,  T2,DT,  AN,  AT,  BN,  BT,  CN,  CT,  IVD,  LAMTAU,  LAMNEB,  CNE 
REAL  TAU (0: 150) , NEB (0:150) ,TDOTl (0:150) ,NDOTl (0:150) , ENE 
REAL  ECHI, CCHI, EALPHA,CALP HA, MU I (0:150) , I, DZ, TCHAR, LCCHI 
REAL  TE, TC, DTAUNDZ, PN, ENR, NR, EGNDB, ELOSSB, RE, SMR, LAMDAR 
REAL  RMUR, KN,NNR, ARECN, PI, CA, CSAHA, NNB (0 : 150) , TAUN (0:150) 

REAL  MSOURCE (0: 150) ,ESOURCE (0:150) ,CV(0:150) ,MUEA(0:150) 

REAL  NDOT 2 (0 : 150) , TDOT2 (0 : 150) , TTILDA (0 : 150) , NTILDA (0 : 150) 

REAL  J,IIVD,AVD,FYEN, EGRADA, PHIB, JNET 


o  o 


INTEGER  N/ CFIX, EFIX, IDEN, CHKDOT, NSTEPS 
C 

COMMON  /PRED/  CA, CSAHA, DT, DTAUNDZ, DZ, EGNDB, ELOSSB, ENR, I, IDEN, KN 
COMMON  /PRED/  LAMDAR, LAMNEB, LAMTAU, MUI, N, NNB, NNR, NR, PI , RE,  RMUR 
COMMON  /PRED/  TAUN, EFIX, CFIX, CHKDOT, FYEN, ENE, ECHI , CCHI , CNE 
COMMON  /PRED/  EALPHA, CALPHA, CV, ESOURCE, LCCHI 

*****HANDLE  EXCEPTIONAL  CONDITIONS 
DZ*1.0/ (N-l) 

DT*(T2-T1) /NSTEPS 
CFIX=0 
EFIX*0 

C*****SET  NEUTRAL  TEMPERATURE  AND  DENSITY 
IF  (TE.EQ.TC) THEN 
DO  10  11*0, N+l 
TAUN (II) *1 . 0 
10  CONTINUE 

ELSE 

DO  20  11*0, N+l 

TAUN (II) *1 . 0+ (TC/TE-1 . 0) *  <11-1.0) / (N-l) 

20  CONTINUE 

END  IF 

NNR*965.5E16*PN/TE 
DO  30  11*0, N+l 

NNB (Il)*1.0/ TAUN (II) 

30  CONTINUE 

DTAUNDZ-TAUN (N) -TAUN  ( 1 ) 

C*****SET  TRANSPORT  PARAMETERS 
RMUR-LAMDAR* SMR 
DO  40  11*0, N+l 

MUI (II) -SQRT (TAUN (II) ) 

40  CONTINUE 

C*****S£T  IONIZATION  AND  SAHA  PARAMETERS 

CA*0 . 41283*ARECN*TCHAR* (NR/1.0E14) **2* (TE/1500) ** (-4.5) 
CSAHA=LOG ( (1 . 4027E20*NNR/NR/NR) *  (TE/1500) **1.5) 


DO  70  ICOUNT*0, NSTEPS-1 

C - PREDICTOR  STEP 

IF  (CHKDOT. EQ. 3)  WRITE (9, 81) 

CALL  DOT (NDOT1, TDOT1, NEB, TAU, EGRADA, PH IB, JNET) 

DO  45  11*0, N+l 

NTILDA (II) *NEB (II) +AN*DT*NDOTl (II) 

TTILDA(Il) *TAU (II) +AT*DT*TDOTl (II) 

45  CONTINUE 

C - CORRECTOR  STEP 

IF  (AN. EQ. 0.0. AND. AT. EQ. 0.0)  THEN 
DO  50  11*0, N+l 
NDOT2 (II) *0 . 0 
TDOT2 (II) *0 . 0 
50  CONTINUE 

GOTO  55 
ELSE 

IF  (CHKDOT. EQ. 3)  WRITE (9, 82) 

CALL  DOT (NDOT2, TDOT2 , NTILDA, TTILDA, EGRADA, PHIB, JNET) 
END  IF 

55  DO  60  11*0, N+l 

NEB (II) *NEB (II) +DT* (BN*NDOTl (II) +CN*NDOT2 (II) ) 

TAU (II) *TAU (II) +DT* (BT*TDOTl (II) +CT*TDOT2  (II) ) 

60  CONTINUE 

70  CONTINUE 


- - 

C*****UPDATE  TIME  DERIV.S,  IMAGE  POINTS,  AND  FIND  PLASMA  POWER  GAIN 
CV(0)=0.0 
CV (N+l)  *0 . 0 
ESOURCE ( 0 ) *0 . 0 
ESOURCE (N+l) *0 . 0 
IF  (CHKDOT.EQ.3)  WRITE (9, 83) 

CALL  DOT (NDOT1 , TDOT1 , NEB, TAU, EGRADA, PHIB, JNET) 

IIVD-0.0 

DO  75  11=2, (N-l) 

IIVD-IIVD+ (ESOURCE (II) -CV(I1) *TD0T1 (II) ) 

75  CONTINUE 

AVD-0.5* (ESOURCE (1) -CV(1) *TDOTl (1) ) +0.5* (ESOURCE (N) -CV(N) 

+  *TDOTl (N) ) 

IVD= (AVD+IIVD) *DZ+2 .0*1* (TAU (1) -TAU (N) ) - (NEB (1) *ENE/KN) 

+  * (TAU (1) -1) 

RETURN 

C 

81  FORMAT (//, ' CALL  DOT  FOR  FIRST  TIME .',//) 

82  FORMAT ( /  / , ' CALL  DOT  FOR  SECOND  TIME.',//) 

83  FORMAT (//, 'CALL  DOT  FOR  LAST  TIME.',//) 

END 


C  DOS  FILE  PRED2.FOR 

C* ************************************************************************** 
SUBROUTINE  DOT (NEBDOT, TAUDOT,  NEB, TAU, EGRADA, PHIB, JNET) 

Q* ********************* ***************************************************** 

REAL  A/ ALPHA (0 : 150) ,BETA(0:150) , CA, CALPHA, CCHI, CDETA, CNE, SIGMA 
REAL  CTETA, CTZ, CV (0 : 150) , CONVECT, D21B, D32B, DT, DZ, DGDU, DELU 
REAL  DETAP, DTAUNDZ, ESOURCE (0:150) , ELOSSB, EGNDB, ETZ , ETETA, DELTAU 
REAL  ENE, ECHI, ENR, F (15) , FYEN, GAMMAP, GAMMAM, GU, I, IB, K  (0:150) 

REAL  LAMTAU, LAMDAR, LAMNEB, MUI (0:150) , MSOURCE (0 : 150) ,MUISOLD 
REAL  MURSOLD, MUIS, MURS, MUEA (0 : 150) , NR, NNR, NEB (0:150) ,NBC0A 
REAL  NBC1A, NBC1B, NBC1C, NES2, NNB (0:150) , NUE, NA (0 : 150) ,NB(0:150) 

REAL  NEBU (0:150) , ND (0: 150) , NS (0 : 150) , NEBDOT (0 : 150) , NEBA (0 : 150) 

REAL  NEBV (0 : 150) , PC (0 : 150) , PI, PB, PBP, POHMIC, P0, QKP, QKM, RE, RMUR 
REAL  TAU (0 : 150) ,TAUN (0:150) , TAUDOT (0 : 150) , TA (0 : 150) , TB  (0 : 150) 

REAL  TS (0 : 150) , TAUA (0: 150) , TAUB (0:150), TAUC (0:150), TAUV(0 : 150) 

REAL  TBC0B, TBC0C, TBC1A, TBC1B, TBC1C, U, CSAHA, DETA, EALPHA, KN, KDE 
REAL  NEBB (0:150) , TAUU (0:150)  ,  NEBC (0:150) , NC (0 : 150) , NBC0B, LCCHI 
REAL  NBC0C, TC (0 : 150) , TBC0A, JNET, SMR, EGRADA, PHIB, BIGU, JRC, CGRADA 
INTEGER  CFIX, EFIX, CHKDOT , CFLAG 
C 

COMMON  /PRED/  CA, CSAHA, DT, DTAUNDZ, DZ, EGNDB, ELOSSB, ENR, I, IDEN, KN 
COMMON  /PRED/  LAMDAR, LAMNEB, LAMTAU, MUI, N, NNB, NNR, NR, PI, RE, RMUR 
COMMON  /PRED/  TAUN, EFIX, CFIX, CHKDOT, FYEN, ENE, ECHI, CCHI, CNE 
COMMON  /PRED/  EALPHA, CALPHA, CV, ESOURCE, LCCHI 
C 

F(l)»5.74E-3 
F  (2) *1 . 40E-3 
F(3)«2.3 
F(4)«0.2 
F  (5) *2 . 70E-2 
F (6) -5 . 74E-3 
F  (7)  -4 . 24E-2 
F (8) *2 . 82 
F (9) *0 . 0 
F (10) *11 . 607 
F (11) *0 . 0 
F (12) *27 .04 

C*****SET  THERMAL  &  ELECTRICAL  CONDUCTIVITIES  AT  0+  (E)  &  1-  (C) 

IF (TAU (1) .LT. 0.1) THEN 
TAU  (1)»0.1 
EFIX-1 
END  IF 

IF  (TAU  (N)  .LT  .0.1)  THEN 
TAU  (N)  -0 . 1 
CFIX-1 
ENDIF 

IF (RE. EQ. 0.5) THEN 
DO  10  I1«0,N+1 
MUEA (II) -TAUN (II) 

10  CONTINUE 

ENDIF 

IF (RE . EQ . 0 . 0 ) THEN 
DO  20  I1-0,N+1 

MUEA (II) "TAUN (II)/ SQRT ( TAU (II) ) 

20  CONTINUE 

ENDIF 

IF (RE. EQ. -0.5) THEN 
DO  30  11-0, N+l 
MUEA (II) -TAUN (II) /TAU(Il) 


30  CONTINUE 

ELSE 

DO  40  11=0/ N+l 

MUEA (II) =TAUN (II) * (TAU(Il) ** (RE-0.5) ) 

40  CONTINUE 

END  IF 

DO  50  11=0, N+l 

K(I1)=( (RE+2.0) /FYEN) *MUEA(I1) *NEB(I1) *TAU(I1) 

PC (II) =NEB (II) * (TAU (II) +TAUN (II) ) 

50  CONTINUE 

DETA=ALOG (K (2) /K (1) ) *DZ/ (K (2) -K (1) ) 

DETAP=ALOG (K (2) /K(l) ) *DZ/ (K (2) -K(l) ) 

C  *  *  *  * ‘DETERMINE  EMITTER  SHEATH***** 

JNET»(I*KN*1 .595769) / (SQRT (TAU (1) ) *NEB  (1) ) 

CALL  SHEATH (JNET, ENR/NEB (1) , TAU (1) , ECHI, PHIB, EALPHA, ENE) 

C - FIND  EMITTER  (0+)  DERIVATIVES  FROM  B.C. 

IF (ECHI . LE . IE-5 . OR . ECHI . GE . 2  0 )  EFIX-1 

ETETA= (TAU (1) -1) *ENE*NEB (1) /KN-I* (ECHI-TAU(l) /2) 

ETZ=ETETA/K ( 1 ) 

EPCZ= (SQRT (PI/8/EALPHA) /LAMDAR/KN) *NEB(1) /MUI (1) -I/MUEA(1) 
ENZ=  (EPCZ-NEB(l)  *  (ETZ+DTAUNDZ)  )  /  (TAU(l)  +TAUN(1)  ) 

C***** SOLVE  COLLECTOR  SHEATH 
CFLAG-0 
CNE=0 . 0 
0-1.0 

GU«G  (U,NEB  (N)  ,  TAU  (N)  ,  I,  KN) 

DELU=0 .05 
DO  80  11=1,50 

DGDU* (G (U+DELU, NEB (N) , TAU (N) , I, KN) -GU) /DELU 

DELTAU=-GU/DGDU 

U=U+DELTAU 

GU=G  (U,  NEB  (N)/TAU(N)  ,I/KN) 

IF  (ABS(GU) .LE. 0.001)  THEN 
CCHI=U*TAU (N) 

GOTO  85 
END  IF 
80  CONTINUE 

CCHI  =0.0 
WRITE (*, 205) 

IF  (CHKDOT  .GT.  1)  WRITE (11, 205) 

C****DETERMINE  DERIVATIVES  AT  COLLECTOR  (1-)  FROM  B.C. 

85  IF  (CHKDOT  .GT.  1) 

4 WRITE (11,206)  EGRADA, ALPHA2 , TAU ( 1 ) , ENR/NEB (1) , 

& JNET, ECHI , PHIB, ENE , EALPHA, CCHI , CALPHA, IFLAG,  RCODE 
DPH  =  -CCHI /TAU (N) 

IF(DPH.LT. 0.0000001)  DPH  »  0.0000001 

CALPHA=(1.0/TAU(N))* (3.14159265/2) *(  (1.0  +  ERF (SQRT (DPH) ) ) 
+  -  (1./2.) * (1 .0  +  ERF (SQRT (4 . 0*DPH) ) )  ) **2  /  (  EXP ( -DPH) 

+  -  ( 1 . / 4 . ) *EXP (-4 . 0*DPH)  ) **2 
LCCHI  -  0.0 
IF (CCHI . LT . 0 . 0 )  THEN 
LCCHI  =  CCHI 
CCHI  =0.0 
ENDIF 

CTETA—I*  (CCHI -TAU  (N)  /2 . 0) 

CTZ-CTETA/K (N) 

CDETA=ALOG(K(N)  /K(N-l)  )  *DZ/  (K  (N)  -K  (N-l)  ) 

NBC0A-TAU ( 1 ) +TAUN ( 1 ) 

NBC0B-SQRT (PI /EALPHA/ 8.0) /LAMDAR/KN /MU I (1) -ENE*NEB (1) /K(l) * 
+  (TAU (1) -1 . 0) -DTAUNDZ 


/aS" 


NBCOOI*NEB(1)  /K (1)  *  (ECHI-TAU(l)  /2.0)  -I/MUEA(1) 

NBC1A=TAU  (N)  +TAUN  (N) 

NBC1B— SQRT  (PI/CALPHA/8 . 0)  /LAMDAR/KN/MUI  (N)  -CNE*NEB  (N)  /K  (N)  * 
+  (TAU  (N)  “1 . 0)  +DTAUNDZ 

NBC1C»-I*NEB (N) /K(N) * (CCHI-TAU(N) /2.0) -I/MUEA(N) 

TBC0A-1.0 

TBC0B»ENE*NEB (1) /KN+I/2 . 0 
TBC0C=*-ENE*NEB  (1)  /KN-I*ECHI 
TBC1A“-1 . 0 

TBC1B=CNE*NEB (N) /KN-I/2 . 0 
TBC1C=-CNE*NEB (N) /KN+I*CCHI 
NEBA  (0) * (NBC0A*LAMNEB) / (2 . 0*DZ) 

NEBB (0) »-LAMNEB*NBC0B 

NEBC (0) *-NBC0A*LAMNEB/ (2 . 0*D2) 

NEBV(0)»NBC0B* (1.0-LAMNEB) *NEB (1) +NBC0C- (1.0-LAMNEB) * 

+  NBCOA* (NEB (2) -NEB (0) ) /  (2.0*DZ) 

TAUA ( 0 ) ®»TBCOA*  LAMTAU / (2.0*DETAP) 

TAUB  (0)  =»-LAMTAU*TBCOB 

TAUC  (0)  — TBCOA*LAMTAU/  (2 . 0*DETAP) 

TAUV(0)»TBC0B*  (1.0-LAMTAU)  *TAU  (1)  +TBCOC-  (1.0-LAMTAU)  * 

+  TBCOA* ( TAU ( 2 ) -TAU (0))/(2.0*DETAP) 

NEBA(N+1) » (NBC1A*LAMNEB) / (2.0*DZ) 

NEBB (N+l ) =-LAMNEB*NBClB 

NEBC  (N+l) «-NBClA*LAMNEB/ (2.0*DZ) 

NEBV  (N+l) *NBC1B* (1.0-LAMNEB) *NEB (N) +NBC1C- (1 . O-LAMNEB) * 

+  NBC1A* (NEB (N+l ) -NEE (N-l ) )/(2.0*DZ) 

TAUA (N+l) -TBC1A*LAMTAU/ (2 . 0*CDETA) 

TAUB (N+l) »-LAMTAU*TBClB 

TAUC  (N+l)  — TBC1A*LAMTAU/  (2 . 0*CDETA) 

TAUV  (N+l) =TBC1B* (1 . 0 -LAMTAU) *TAU (N) +TBC1C- (1 . Q-LAMTAU) * 

+  TBC1A* (TAU (N+l) -TAU (N-l) ) / (2 . 0*CDETA) 

C***** INITIALIZE  GAMMAP  &  QKP  FOR  LOOP 

MURS-MUI (2) /MUEA (2) +MUI (1) /MUEA(l) *  <1-2*DZ*  ( 

+  (0.5-RE) *ETZ/TAU (1) -0 ,5*DTAUNDZ/TAUN (1) ) ) 

GAMMAP-0.5*  ( ( (MUI  (1)+MUI  (2)  )  *  (PC(1)  -PC(0) ) 

+  )/DZ+I*MURS) 

QKP®*  (TAU  ( 1 )  -TAU  ( 0 ) )  /DETA 
MUIS=MUI (1) +MUI (2) 


DO  100  J*1#N 

C - UPDATE  FOR  NEW  J 

GAMMAM-GAMMAP 

QKM-QKP 

DETA*DETAP 

MUISOLD'-MUIS 

MURSOLD-MURS 

IF  (J.NE.N)  THEN 

DETAP-ALOG(K(J+l)  /K(J)  )  *DZ/  (K(J+1)-K(J) ) 

MUIS-MUI (J) +MUI (J+l) 

MURS-MUI (J) /MUEA ( J) +MUI (J+l) /MUEA (J+l) 

ELSE 

MURS-MURS+2+DZ* (MUI (N) /MUEA(N) ) * ( (0.5-RE) *CTZ/TAU(N) 
+  -0.5*DTAUNDZ/TAUN(N) ) 

END  IF 

C - FIND  AMBIPOLAR  FLUX  AT  J+l/ 2 

GAMMAP-0.5* ( (MUIS* (PC (J+l) -PC (J) ) ) /DZ+I*MURS) 

C - FIND  MASS  SOURCE  AT  J 

A-CA/TAU(J)  **4.5 

NES2-NNB(J)  *TAU  ( J)  **1 ,5*EXP  (CSAHA-EGNDB/TAU  (J)  ) 

D21B-F (7) * (1+F (8) /TAU ( J) ) 


D32B-F (2) *EXP (F (3) /TAU  ( J)  ) 

>  IB=A*NES2* (1+F (1) /NEB ( J) ) / (1+D21B* (1+D32B/NEB  ( J)  ) /NEB  ( J)  ) 

P0-1+  (F  (4)  /NEB  (J)  )  *  (1+F  (5)  /NEB  (J)  )  /  (1+F  (6)  /NEB  (J)  ) 

NUE-NEB (J) *NEB (J) /NES2 
MSOURCE ( J) -NEB (J) *IB* (1-P0*NUE) 

IF  (IDEN.EQ.l)  MSOURCE (J) -NEB (J) *A*NES2 
NEBDOT  ( J)  -RMUR*  (GAMMAP-GAMMAM)  /DZ+MSOURCE  ( J) 
NA(J)-RMUR*MUIS*  (TAU(J+1)  +TAUN(J+1)  )  /2.0/DZ**2 
NB (J) -RMUR* (MUIS+MUISOLD) * (TAU ( J) +TAUN ( J) ) /2 . 0/DZ**2 
NC ( J) =RMUR*MUISOLD* (TAU (J-l) +TAUN (J-l) ) /2.0/DZ**2 
ND (J) -I* (MURS-MURSOLD) *RMUR/DZ/2 . 0 
+  +NEB ( J) *IB* (1.0+SQRT (P0/NES2) *NEB(J) ) 

C  NS ( J) -IB* (1 . 0-PQ*NUE) 

NS  ( J)  —NEB  ( J)  *IB*  (1 . 0+SQRT  (P0/NES2)  *NEB  ( J)  )  *SQRT  (P0/NES2) 

IF  (IDEN.EQ.l)  NS ( J) =A*NES2 
NEBA(J)— DT*NA(J)  *LAMNEB 

NEBB ( J) -1 . 0+DT*NB ( J) *LAMNEB-DT*NS ( J) *LAMNEB 
NEBC ( J) — DT*NC ( J) *LAMNEB 

NEBV  ( J) —NEB ( J) +DT*NA ( J) * (1. 0-LAMNEB) *NEB ( J+l) -DT*NB ( J) * 

+  (1. 0-LAMNEB) *NEB (J) +DT*NC (J) * (1 . 0-LAMNEB) *NEB (J-l) + 

+  DT*ND (J) +DT*NS (J) * (1 . 0-LAMNEB) *NEB ( J) 

KDE-K ( J) * (DETA+DETAP) / 2 
QKP- (TAU (J+l) -TAU ( J) ) /DETAP 

CONVECT—  (1.5)  *1*  ( DETA*  QKP + DETAP  *  QKM )  /  (2*KDE) 

SIGMA-NEB  ( J)  *MUEA(  J) 

POHMIC-I*  (I/SIGMA+TAU  (J)  *  (NEB  (J+l)  -NEB  (J-l)  ) 

+  / <2*DZ*NEB(J) ) ) 

PBP— (F (9) *NNR/NR) *EXP (-F (10) /TAU ( J) ) 

PB- (F (11) *NNR/NR) *EXP (-F (12) /TAU ( J) ) 

CV (J)  -1 .5*NEB (J) +NNB ( J) * (F(10) *PBP+F (12) *PB*NUE) 

+  /  (TAU  ( J)  *TAU  ( J)  ) 

ESOURCE ( J) — ELOSSB*MSOURCE ( J) 

+  -NNB ( J) *PB* (2*NUE*NEBDOT ( J) /NEB (J) ) 

TAUDOT ( J) « ( (QKP-QKM) /KDE+CONVECT+POHMIC+ESOURCE ( J) ) 

+  /CV(J) 

TA( J) -1 .0/ (DETAP *KDE*CV (J) ) 

TB ( J) » (1 . 0/DETAP+1 . 0/DETA) /KDE/CV  ( J) 

TC  ( J) -1 . 0/DETA/KDE/CV ( J) 

TS ( J) - (CONVECT+POHMIC+ESOURCE ( J) ) /CV ( J) 

TAUA(J) — DT*LAMTAU*TA ( J) 

TAUB ( J) -1 . 0+DT*LAMTAU*TB ( J) 

TAUC  ( J)  — DT*LAMTAU*TC  ( J) 

TAUV(J)  -TAU  ( J)  +DT*  (1.0-LAMTAU)  *TA(J)  *TAU(J+1) 

+  -DT* (1.0-LAMTAU) *TB ( J) *TAU ( J) 

+  +DT* (1.0-LAMTAU) *TC ( J) *TAU ( J-l) 

+  +TS (J) *DT 

IF  (CHKDOT.EQ.3)  WRITE (9, 201)  J, NEB ( J) , J, TAU ( J) , 

+  J/ MSOURCE (J) , J,PB,PBP,A,D21B,D32B,P0, IB, 

+  NUE,NES2, QKP, GAMMAP, DETAP, MURS 

100  CONTINUE 

- - 

NEBC (0) -NEBC (0) -NEBC (1) *NEBA (0) /NEBA(l) 

NEBB (0) -NEBB (0) -NEBB (1) *NEBA (0) /NEBA(l) 

NEBV ( 0 ) -NEBV ( 0 ) -NEBV ( 1 ) *NEBA ( 0 ) /NEBA ( 1 ) 

NEBA (0) -NEBB (0) 

NEBB (0) -NEBC (0) 

NEBA(N+1) -NEBA (N+l) -NEBA (N) *NEBC (N+l) /NEBC  (N) 

NEBB (N+l) -NEBB (N+l) -NEBB (N) *NEBC  (N+l) /NEBC  (N) 

NEBV  (N+l) -NEBV (N+l) -NEBV (N) *NEBC(N+1) /NEBC (N) 

NEBC (N+l) -NEBB (N+l) 
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NEBB (N+l) =NEBA (N+l ) 

TAUC ( 0 ) -TAUC { 0 ) -TAUC ( 1 ) *  TAUA ( 0 ) / TAUA ( 1 ) 

TAUB (0) -TAUB (0) -TAUB (1) *TAUA(0) /TAUA(l) 

TAUV ( 0 )  -TAUV ( 0 )  -TAUV ( 1 )  *TAUA(0)  /TAUA(l) 

TAUA ( 0 ) -TAUB ( 0 ) 

TAUB ( 0 ) -TAUC ( 0 ) 

TAUA  (N+l)  -TAUA (N+l)  -TAUA  (N)  *TAUC  (N+l)  /TAUC  (N) 

TAUB  (N+l)  -TAUB  (N+l)  -TAUB  (N)  *TAUC  (N+l)  /TAUC  (N) 

TAUV  (N+l)  -TAUV  (N+l)  -TAUV(N)  *TAUC  (N+l)  /TAUC  (N) 

TAUC (N+l) -TAUB (N+l) 

TAUB (N+l ) -TAUA (N+l ) 

ALPHA (0)— NEBA(O)  /NEBB  (0) 

DO  110  11-1/ N 

ALPHA (II) — NEBA (II)/ (NEBC (II) * ALPHA (11-1) +NEBB  (II) ) 
110  CONTINUE 

ALPHA (N+l) -0.0 

BETA ( 0 ) -NEBV ( 0 ) /NEBB  ( 0 ) 

DO  120  Il-l/N+1 

BETA  (ID  —  (NEBV  (II)  -NEBC  (II)  *BETA(I1-1)  )  /  (NEBC  (II) 

+  *ALPHA(I1-1) +NEBB (II) ) 

120  CONTINUE 

NEBU (N+l) -BETA (N+l) 

DO  130  Il-N,0/-1 

NEBU  (ID -ALPHA  (II)  *NEBU  (1 1  +  1)  +BETA(I1) 

130  CONTINUE 

ALPHA  (0)  --TAUA(O)  /TAUB(0) 

DO  140  Il-l/N 

ALPHA  (ID—  TAUA  (II)/  (TAUC  (II)  *  ALPHA  (11-1)  +TAUB  (II)  ) 
140  CONTINUE 

ALPHA (N+l) -0.0 

BE  T A ( 0 ) - TAUV ( 0 ) / TAUB ( 0 ) 

DO  150  11-1,  (N+l) 

BETA  (ID  —  (TAUV  (II)  -TAUC  (II)  *BETA  (11-1) )/  (TAUC  (II) 

+  *ALPHA(I1-1) +TAUB (II) ) 

150  CONTINUE 

TAUU (N+l) -BETA (N+ 1 ) 

DO  160  Il-N/0,-1 

TAUU  (ID -ALPHA  (II)  *TAUU  (11+1) +BETA  (II) 

160  CONTINUE 

DO  170  II— 0, N+l 

TAUDOT (II) » (TAUU (II) -TAU(Il) ) /DT 

NEBDOT  (II)  *  (NEBU  (ID  -NEB  (II) )  /DT 
170  CONTINUE 

IF  (CHKDOT.EQ. 3)  THEN 

WRITE (9/ 202)  (II,  NEBDOT (II) ,11, TAUDOT (II) ,11-0/ N+l) 

WRITE (9,203)  (II, ALPHA (II) , II, BETA (II) , I1-0,N+1) 

WRITE (9,204)  (II, NEBU (II) ,11, TAUU (II) ,11-0, N+l) 

END  IF 
RETURN 

* 

201  F0RMAT('NEB(',I2,')-',F8.3,'  TAU  (',12, ' )  »' ,  F8 . 3, 

6'  MSOURCE(',I2,')-',F8.3/'  J=',I2,  '  PB-',F8.3, 

6'  PBP-',F8.3/'A-',F8.3,'  D21B-',F8.3,'  D32B-',F8.3/ 

6'P0-' ,F8.3, '  IB-' ,F8.3/'NUE»',F8.3,'  NES2-/,F8.3/ 

4'QKP-',F8.3, '  GAMMAP-' , F8 . 3/' DETAP-' , F8 . 3, '  MURS-',F8.3) 

202  FORMAT (' NEBDOT (', 12, ') »'  ,F8. 3, '  TAUDOT (' , 12, ' ) -' , F8 . 3) 

203  FORMAT  (' ALPHA  (' , 12, ' ) -' , F8 . 3, '  BETA  (',12, ')  -' ,  F8 . 3) 

204  FORMAT  ('NEBU  (',  12, ')«',F8. 3,'  TAUU  (',12, ' )  ,  F8 . 3) 

205  FORMAT (//2X, 'COLLECTOR  SHEATH  FAILED  TO  CONVERGE',//) 

206  FORMAT  (/IX, '  EGRADA-' ,  F7 . 3,  3X, '  ALPHA2-' ,  F7 . 3/ 


41X, 'RLE-' , F5 . 2, 3X, ' EMISS-' , F9 . 1 , 3X, ' JNET-' , F7 . 4/ 

41X, » ECHI-' , F7 . 3, 3X, ' PHIB-' , F7 . 3, IX, ' ENE-' , F10 . 3, 3X, 

4' EALPHA-' , F7 . 4 / IX, '  CCHI-' , F7 . 3 , 3X, ' CALPHA-' , F7 . 4 / IX, 

4' IFLAG-' ,11, 3X, ' RCODE-' , A1 ) 

END 

Q****** **********************  **********  ************************************ 

FUNCTION  G  (UX, NEB, TAU, I, KN) 

C********* ***************************************************************** 

REAL  UX,NEB, TAU, I, KN, SQX, GX 
DOUBLE  PRECISION  ERF 
IF (UX.LE.O .0) THEN 
SQX-0 . 0 
ELSE 

SQX-SQRT <UX) 

END  IF 

C  G=UX*LOG (1 .0+ERF (SQX) )  -  LOG (NEB*SQRT (TAU) /I/KN/2 . 0) 

G-UX  -  LOG(NEB*SQRT (TAU) /I/KN/2. 0) 

RETURN 

END 

Q******** ********************************************************** 

FUNCTION  ERF (X) 

Q***************** ******  ******  ************************  ************* 

DOUBLE  PRECISION  XI,  SUM,  ERF,  ERFX, P (3, 6) , Q (3, 6) , PI, TSUM1, TSUM2 
INTEGER  IPOWER(3) 

DATA  IPOWER/2,1,-2/ 

PI-3.1415926535898 
TSUM1-0 . 0 
TSUM2-0.0 
Xl-X 

IF  (XI  .LT.  0.0)  XI— XI 
IF  (XI  .GT.  5.93)  THEN 
ERFX-1 . 0 
GOTO  1240 
END  IF 

IF  (XI  .EQ.  0.0)  THEN 
ERF-0 . 0 
GOTO  1245 
ELSE 

IFLAG-1 
END  IF 

IF  (XI  .LE.  4.0  .AND.  XI  .GE.  0.47)  IFLAG-2 
IF  (XI  .GT.  4.0)  IFLAG-3 
DO  1205  J-1,6 

TSUM1-TSUM1  +  P (IFLAG, J)  *  (XI** (IPOWER(IFLAG) * ( J-l) ) ) 

TSUM2-TSUM2  +  Q (IFLAG, J)  *  (XI** (I POWER (IFLAG) * (J-l) ) ) 

1205  CONTINUE 

SUM-TSUM1 /TSUM2 
GOTO  (1210,1220,1230) , IFLAG 
1210  ERFX-X1*SUM 
GOTO  1240 

1220  ERFX-1. 0  -  DEXP(-X1*X1)*SUM 
GOTO  1240 

1230  ERFX-1. 0  -  DEXP (-X1*X1) /XI* (1 . 0/DSQRT (PI)  +  (1 . 0  /  (XI* XI) *SUM) ) 

1240  ERF -ERFX 

IF  (X  .LT.  0.0)  ERF— ERFX 
1245  RETURN 

q****** **********0 (IFLAG, J) ******Q (IFLAG, J) ****** 

DATA  P (1, 1) /3 .20937758913847D+03/ , P (1, 2) /3.774852376853D+02/ 

DATA  P(l, 3) /1.1386415415105D+02/,P(l, 4) /3. 16112374387057/ 

DATA  P (1,5) / 1.8577770 61 8 4 603D-01/,P (1, 6) /0.0/ 


DATA  Q (1/ 1) /2 . 84423683343917D+03/ / Q (1, 2) /I .2826165077372D+03/ 
DATA  Q(l,3) /2 . 44024 637 9344 44D+02/, Q (1,  4) /2 . 36012909523441D+01/ 
DATA  Q(l,5)/1.0/,Q(l,6)/0.0/,P(2,l)  /2 . 2898992851659D+01/ 

DATA  P(2,2) /2 . 6094746956075D+01/, P (2,3) /I . 4571898596926D+01/ 

DATA  P (2, 4) /4. 2677201070898/, P (2,5) /5 . 6437160686381D-01/ 

DATA  P (2, 6) /-6. 0858151959688D-06/, Q (2 , 1) /2 . 289898574 98 91D+01/ 
DATA  Q (2, 2) /5 . 1933570 687 552D+01/, Q (2,3) /5 . 02732028 63803D+01/ 

DATA  Q (2, 4) /2 . 6288795758761D+01/, Q (2,5) /7 .5688482293618/ 

DATA  Q (2, 6) /1.0/,P (3,1) /-6 . 5874916152 9838D-04/ 

DATA  P (3, 2)  /-l. 6083785148742 3D-02/ ,P (3,  3) /-1.2578172611123D-01/ 
DATA  P (3,  4) /-3.60344899949804D-01/,P (3,5) /-3 . 05326634961232D-01/ 
DATA  P (3, 6)/-l. 6315387137302 lD-02/,Q(3, 1) /2 . 3352049762687D-03/ 
DATA  Q (3, 2) /6 . 05183413124413D-02/ , Q (3,3) /5 .27905102951428D-01/ 
DATA  Q (3, 4) /1 . 87295284992346/, Q (3, 5) /2 . 56852019228 982/, Q  (3, 6) / 1 / 
END 


CALPHA-0 . 5 
DO  10  I1»0,N+1 
*  NEB (II) =11 

NEBCAL (II)* (NEB (II) -1) / (N-l) 

NEB (II) “4 . 0* (KN+NEBCAL (II) * (1-NEBCAL (II) ) ) 
TAU(Il)-2700/TE  • 

10  CONTINUE 
END  IF 

c*****miscelaneous  DEFAULTS 
ARECN-0.31 

EGNDB-3. 896/8. 609E-05/TE 
ELOSSB-EGNDB 

SN-DELTAT* (N-l) **2*3*LAMDAR*SMR*  (BN+CN) 

ST-DELTAT* (N-l) **2*0.667* (RE+2) *2** (RE+0.5) * (BT+CT) /FYEN 
C*****NONDIMENSIONALIZE  CURRENT  AND  CALCULATE  RICHARDSON  EMISSION 
VALUE-EXP (-11600 . 0*EWF/TE) 

ENR- (7 . 676E+14* (TE) **1 . 5*VALUE) /NR 
I-J  /  (KN*NR* (3 . 1265322E-13) *SQRT (TE) ) 

JRIC-120 . 0*TE*TE* (EXP (-11600 . 0*EWF/TE) ) 

C*****0UTPUT  INITIALIZATION  DATA  TO  PRNTOUT.DAT 

WRITE (9/ 155)  TE, TC,EWF, CWF, PN, D, J, CHKDOT, OFILE,N, 

&  JRIC, NR, TCHAR, I , ENR, KN, SMR, LAMDAR, 

5  NSTEPS,T2, DELTAT, DTP, LS 
IF  (CHKDOT. GT.l)  THEN 

WRITE (9, 160)  ECHI,ENE,EALPHA, CCHI, CNE, CALPHA, AN, AT, 

6  BN,  BT,  CN, CT, SN, ST, LAMNEB, LAMTAU, ELOSSB, 

&  ARECN, EGNDB, IDEN, FYEN, RE 

WRITE (9, 159)  (II, NEB (II), II, TAU (II) , II— 1,N) 

END  IF 

WRITE (9,161) 

C*****CHECK  FOR  UNREAL  CURRENT  CONDITION  J/JR  >1.0 
IF  (2*KN*I/ENR  .GE.  1.0)  THEN 
WRITE (*,153) 

STOP 

ENDIF 

RETURN 


150  FORMAT (IX, '*****************************************'  / 

SIX, '  *****  *★***'/ 

SIX, ******  TEC  START  *****'/ 

SIX,'*****  *****'/ 

SIX,'*****************************************'///) 

151  FORMAT (F8 i 1/F8.1/F6.3/F6 . 3/F6 . 3/F6 . 3/F7 . 2/F5 .2/F6 . 1/ 
SI1/I1/I3/I3/I3) 

152  FORMAT (F9.4/F9.4/F8.3/F8.3/F9.4/F9. 4/13) 

153  FORMAT (IX, '***SMALL  J  IS  TOO  LARGE... CASE  TERMINATED***'//) 

154  FORMAT (2F8. 3) 

155  FORMAT (12X, '  TEC  INITIAL  DATA  SUMMARY'/ 


SIX, 

SIX, 

six, 

SIX, 

SIX, 

SIX, 

six, 

six, 

six, 

six, 

SIX, 


PHYSICAL  OPERATING  CONDITIONS 
EMITTER  TEMPERATURE 
COLLECTOR  TEMPERATURE 
EMITTER  WORK  FUNCTION 
COLLECTOR  WORK  FUNCTION 
CONVERTOR  PRESSURE 
GAP  THICKNESS 
OPERATING  CURRENT 
TEC  FUNCTION  SETTINGS 


- -  n 

— '// 

(TE) -' ,F8.1, '  KELVIN'/ 

(TC) ■' , F8 . 1, '  KELVIN'/ 

(EWF) -  ' , F6 . 3, '  EV'/ 

(CWF)-  '  ,  F6 . 3, '  EV'/ 

(PN)  »  ' ,  F6 . 3, '  TORR'/ 

(D) —  ' ,F6.3, '  MM'/ 

(J)-  '  ,  F7 . 3, '  AMP  S  /  CMA  2 '  /  / 
. '// 


DIAGNOSTIC  LEVEL  (CHKDOT) 1 


,11/ 


c  DOS  FILE  PEED3 .FOR 

£♦********************************************************************* 
SUBROUTINE  INITIAL (TE, TC, EWF, CWF, PN, NSTEPS, DTP, T2, AN, AT, 

&CN, CT,BN, BT, TCHAR, SMR/ ARECN/ DELTAT, SN, ST, TAU, NEB, LS, pc) 

C* ********************************************************************* 

REAL  CA, CSAHA, CNE, ENE, ECHI, CCHI, EALPHA, CALPHA, LAMNEB, LAMTAU 
REAL  DT, DTAUNDZ, DTP, T2, AN, AT, CN, CT, BN, BT, RE, KN, TCHAR, PN 
REAL  SMR, LAMDAR, NR,  TE, TC, ENR, I , ARECN, EGNDB, ELOSSB,MUI (0:150) 

REAL  TAU (0:150) ,NEB(0:150) ,DELTAT, SN, ST,PI,TAUN (0 :150) ,NNB(0:150) 
REAL  NEBCAL  (0:1-50)  ,  FYEN ,  LAMDAE ,  LAMDAI ,  EWF ,  CWF ,  J ,  RMUR,  LCCHI 
REAL  ESOURCE (0:150) , CV, JRIC, VALUE, PHISDAT (154, 6) , PHIBDAT (21, 6) 
INTEGER  N, IDEN, CHKDOT, OFILE, EFIX, CFIX, NSTEPS, LS, pc 
C 

COMMON  /PRED/  CA, CSAHA, DT, DTAUNDZ, DZ, EGNDB, ELOS SB, ENR, I, IDEN, KN 
COMMON  /PRED/  LAMDAR, LAMNEB, LAMTAU, MUI, N, NNB, NNR, NR, PI, RE, RMUR 
COMMON  /PRED/  TAUN, EFIX, CFIX, CHKDOT, FYEN, ENE, ECHI, CCHI, CNE 
COMMON  /PRED/  EALPHA, CALPHA, CV, ESOURCE, LCCHI 
COMMON  /XSHEATH/  PHISDAT, PHIBDAT 
C 

WRITE (*, 150) 

C*****reAD  FILE  INDATA . DAT 

READ (2, 151)  TE, TC, EWF, CWF, PN,D,J, DTP, T2, CHKDOT, OFILE, pc, N,LS 
REWIND (2) 

CLOSE (2) 

C*****READ  FILE  PRECOR.DAT 
C  CALL  DATAINT 

C*****SET  NUMERICAL  PARAMETERS,  RECOMPILATION  REQUIRED  TO  CHANGE 
AN-0.5 
AT-0.5 
CN-0.5 
CT-0.5 
BN-1.0-CN 
BT-1.0-CN 
PI-3.1415926 
IDEN-0 
RE-0.0 
FYEN-1 . 0 
NR-1.0E14 
DELTAT-DTP /NSTEPS 
LAMNEB- 1 . 0 
LAMTAU- 1 . 0 

C*****SET  TRANSPORT  PROPERTIES 
LAMDAE-1 . 0/32 . 3/PN 
LAMDAI-1 . 0/96 . 6/PN 
LAMDAR- LAMDAI / LAMDAE 
KN-LAMDAE/D 

TCHAR-D/ (KN*3 . 75* (TE**.5) ) 

SMR-1.0/492.2 

C*  *  *  * ‘EXECUTE  OFILE  SELECTION  SETTING 
IF  (OFILE  .EQ.  1)  THEN 

READ (8,152)  ENE , CNE , ECHI , CCHI , EALPHA, CALPHA, N 
READ (8, 154)  (NEB (II) , TAU (II) , I1-0,N+1) 

REWIND (8) 

ELSE 

ENE-0 . 8 
CNE-0 . 8 
ECHI-3.0 
CCHI-3.0 
EALPHA-0 . 5 


)2>l 


SIX,'  RESTART  SEQUENCE  (OFILE) *  ',11/ 

SIX,'  POINT  DENSITY  (N)*',I3// 

SIX, 'PHYSICAL  PARAMETERS  EVALUATED - '// 

SIX,'  RICHARDSON  CURRENT  ( JRIC) *' , E9 . 2, '  AMPS/CM*2' / 
SIX,'  REFERENCE  DENSITY  (NR) *' ,E9.2, '  1/CMA3' / 

SIX,'  CHARACTERISTIC  TIME  (TCHAR)  «  ',F7.4,'  SECS*E-06'/ 

SIX,'  NONDIM  CURRENT  (I)*  ',F7.4/ 

SIX,'  NONDIM  EMISSION  (ENR) -  ' ,F8 . 3, ' (NRIC/NR) ' / 

SIX,'  KNUD SEN  NUMBER  (KN) *  ',F7.4/ 

SIX,'  SQRT (MASS  RATIO)  <SMR)  =  ',F7.4/ 

SIX, 'MEAN  FREE  PATH  RATIO  (LAMDAR) -  ',F7.4// 

SIX, 'TIME  SETTINGS - ' /4X, ' NSTEPS*' , 13/ 

S4X, '  T2*',F6.1/4X,'DELTAT=',F6.3/4X,'  DTP-',F6.3/ 

S4X, '  LSF«',I3) 

159  FORMAT  (4X,  'NEB  (' ,  12,  ' ) *' , F8 . 3, '  TAU  (' ,  12, ' )  ='  ,F8 . 3) 

160  FORMAT  (//IX, 'ADVANCED  DIAGNOSTIC  OUTPUT - '// 

S4X, '  ECHI*' ,  F5 . 1, '  ENE*'  ,  F7 . 4,  '  EALPHA=' ,  F7 . 4/ 

S4X,  'CCHI=*' ,  F5.1,  '  CNE*'  ,F7 . 4,  '  CALPHA-' , F7 . 4/ 

S4X, 'AN*' , F7 . 4, '  AT*' , F7 . 4/4X, ' BN-' , F7 . 4, '  BT*',F7.4/ 

S4X, 'CN*',F7.4, '  CT*',F7.4/4X,'SN=',F7.4,'  ST-',F7.4/ 

S4X, 'LAMNEB*  ' ,F5.2/4X, 'LAMTAU*' ,F5.2/4X, 'ELOSSB*  ',F6.3/ 
S4X, ' ARECN*  ' , F6 . 3/4X, ' EGNDB*  ' , F6 . 3, /4X, ' IDEN*' , 11/ 

S4X, ' FYEN*  ',F5.2/4X, 'RE*  ',F5.2// 

SIX, 'STARTUP  DENSITY  AND  TEMPERATURE  RATIOS - ') 

161  FORMAT  (/IX, ' - ', 

- /) 

END 
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C  DOS  FILE  PRED4 .FOR 

SUBROUTINE  SHEATH (JNET, ENR, TAU, ECHI , PHIB, EALPHA, ENE) 

REAL  ENE , ECHI , ENR/ TAU , PHIB , EALPHA, JNET , A, S / B , DPH 
A=1 . 1485 
B=1 .169 
S-.l 

ECHI-( ( (JNET-A) +SQRT ( (JNET-A) **2+4*S*B) ) /  (2*S) ) **2 

IF  (ECHI. LT. 0.0)  ECHI  *  0.0 

ENE=(SQRT (3. 14159265/2) * JNET* SQRT (TAU) )  + 

+  (1/EXP (ECHI/TAU)) 

PHIB=ALOG (ABS (ENR/ENE) ) 

IF (PHIB. LE. 0.0) THEN 
PHIB  -  0.0 
ENE  *  ENR 

ECHI  »  -  TAU*ALOG (ENE  -  (SQRT (3 . 14159265/2) *JNET*SQRT (TAU) ) ) 
END  IF 

DPH  *  (PHIB  -  ECHI) /TAU 

IF(DPH.LT. 0.0000001)  DPH  -  0.0000001 

EALPHA* (1 . 0/TAU) * (3.14159265/2) * (  (1.0  +  ERF (SQRT (DPH) ) ) 

+  -  (1 . /2 . ) *  (1 . 0  +  ERF (SQRT (4 . 0*DPH) ) )  ) **2  /  (  EXP(-DPH) 

+  -  ( 1 . / 4 . ) *EXP (-4.0 *DPH)  ) **2 
RETURN 
END 
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